%% Task 4: Simulated EMA
L=1; width=0.020; height=0.010; rho=7870; E=210e9; Ixx=(width*height^3)/12;
A=width*height; V=L*A; m=V*rho; N=5; xin=1; xout=0:.1:1;
zeta=1; beta=fzero('cos(x).*cosh(x)+1',[1,2]);
f0=sqrt(beta^4.*E.*Ixx./(rho*A))/(2*pi);

fs=10e3;fc=1e3;df=1/fs;
force=randn(512*1024,1);

[B,A]=butter(4,fc/(fs/2));
filtForce=filter(B,A,force);

psdForce=pwelch(filtForce);

figure()
plot(psdForce)
xlim([0 fc])

rms(psdForce)
rms(filtForce)

fftSize=2^14;
for i=1:11
    [poles(:,i),residues(:,i)]=cantilever(m,f0,xin,xout(i),N,zeta);
    [y(:,i),t(:,i)] = timeresv_rp(filtForce,fs,residues(:,i),poles(:,i),(1:N));
    FRF(:,i)=tfestimate(t(:,i),y(:,i),hann(fftSize),(fftSize/2),fftSize);
    figure()
    plot(20*log10(abs(FRF(:,i))));hold on;xlim([0 fc]);
end
hold off
cohere=mscohere(filtForce,y(:,3));
figure()
plot(cohere);xlim([0 fc])

% MIF=modeind1(FRF);
% figure()
% plot(MIF);xlim([0 fc])