modal_frf.zip_modalfrf_transfer

disp(' ') disp(' modal_frf.m ver 1.3 July 20, 2013 ') disp(' by Tom Irvine Email: tom@vibrationdata.com ') disp(' ') disp(' This program calculates the frequency response function ') disp(' and coherence of a force and response time history pair. ') disp(' It is best suited for a random vibration force input.') disp(' It processes the data as ensemble average '); % close all; % clear n; clear Y; clear amp; clear tim; clear N; clear df; clear dt; clear Mag; clear mag_seg; clear amp_seg; clear full; clear freq; clear ss; clear seg; clear i_seg; clear rlab; clear ylab; % clear t1; clear t2; clear A; clear B; % clear PSD_force; clear PSD_resp; % clear H1; clear H1_mag; clear H1_phase; clear H1_m; clear H1_p; % clear H2; clear H2_mag; clear H2_phase; clear H2_m; clear H2_p; % clear COH; % fig_num=1; % disp(' '); disp(' Select response amplitude type: '); disp(' 1=displacement '); disp(' 2=velocity '); disp(' 3=acceleration '); iamp=input(' '); % if(iamp==1) tlabel='Receptance'; end if(iamp==2) tlabel='Mobility'; end if(iamp==3) tlabel='Acceleration'; end % disp(' '); disp(' Select input method: '); disp(' 1=Signals are in two separate arrays. Each with time(sec) & amp '); disp(' 2=Signals are in a common array. time(sec) & force & response '); im=input(' '); % iflag=0; % if(im==1) disp(' '); disp(' Force '); disp(' '); [t1,A,dt1,sr1,tmx1,tmi1,n1,~]=enter_time_history(); disp(' '); disp(' Response '); disp(' '); [t2,B,dt2,sr2,tmx2,tmi2,n2,ncontinue]=enter_time_history(); % if(tmx1 ~= tmx2 || tmi1 ~= tmi2 || n1 ~= n2 ) disp(' Error: time columns must be the same '); iflag=1; else dt=dt1; sr=sr1; t=t1; n=n1; tmx=tmx1; tmi=tmi1; end % else [t,A,B,dt,sr,tmx,tmi,n,ncontinue]=enter_time_history_three(); end % if(ncontinue==1) % disp(' ') disp(' mean removal? '); mr_choice=input(' 1=yes 2=no ' ); disp(' ') % disp(' ') disp(' Select Window ') h_choice = input(' 1=Rectangular 2=Hanning '); % %%%%%%%%%%%%%%% advise % [df,mmm,NW]=PSD_advise(dt,n); % %%% begin overlap % mH=((mmm/2)-1); % full=zeros(mH,1); mag_seg=zeros(mH,1); % nov=0; % clear amp_seg_force; clear amp_seg_resp; % H1=zeros(mmm,1); H2=zeros(mmm,1); % for ijk=1:(2*NW-1) % amp_seg_force=zeros(mmm,1); amp_seg_force(1:mmm)=A((1+ nov):(mmm+ nov)); % amp_seg_resp=zeros(mmm,1); amp_seg_resp(1:mmm)=B((1+ nov):(mmm+ nov)); % nov=nov+fix(mmm/2); % [complex_FFT_force]=CFFT_core(amp_seg_force,mmm,mH,mr_choice,h_choice); [complex_FFT_resp]=CFFT_core(amp_seg_resp,mmm,mH,mr_choice,h_choice); % dH1=(conj(complex_FFT_force)).*complex_FFT_force; % two-sided dH2=(conj(complex_FFT_resp)).*complex_FFT_force; % H1=H1+((conj(complex_FFT_force)).*complex_FFT_resp)./dH1; H2=H2+((conj(complex_FFT_resp)).*complex_FFT_resp)./dH2; % end % den=df*(2*NW-1); H1=H1/den; H2=H2/den; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Impulse Response Function % clear Y; clear amp_real; clear amp_imag; Y = ifft(H1,mmm,'symmetric'); if(iamp==3) Y(1)=0; end amp_real=real(Y); amp_imag=imag(Y); % t=fix_size(t); tt=t(1:mmm); t_real=fix_size(amp_real); t_imag=fix_size(amp_imag); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % H1_mag=abs(H1); H1_phase=(180/pi)*atan2(imag(H1),real(H1)); % H2_mag=abs(H1); H2_phase=(180/pi)*atan2(imag(H2),real(H2)); % H1_m=H1_mag(1:mH); H1_p=H1_phase(1:mH); % H2_m=H2_mag(1:mH); H2_p=H2_phase(1:mH); % clear H1C; clear H2C; % H1C=H1(1:mH); H2C=H2(1:mH); % COH=zeros(mH,1); for i=1:mH COH(i)=abs(H1(i)/H2(i)); end % fmax=(mH-1)*df; freq=linspace(0,fmax,mH); % disp(' '); f1=input(' Enter plot lower frequency(Hz) '); f2=input(' Enter plot upper frequency(Hz) '); % figure(fig_num); fig_num=fig_num+1; subplot(3,1,1); plot(freq,H1_p); out1=sprintf('H1 %s Frequency Response Function',tlabel); title(out1); grid on; ylabel('Phase (deg)'); axis([f1,f2,-180,180]); set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','lin','ytick',[-180,-90,0,90,180]); % subplot(3,1,[2 3]); plot(freq,H1_m); grid on; xlabel('Frequency(Hz)'); ylabel('G_F_X/G_F_F'); set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log',... 'YScale','log'); % ymax=10^(ceil(log10(max(H1_m)))); ymin=ymax; for i=1:mH if(H1_m(i)<ymin && freq(i)>f1 && freq(i)<f2 ) ymin=H1_m(i); end end % if(ymin<ymax/10000) ymin=ymax/10000; end ymin=10^(floor(log10(ymin))); % axis([f1,f2,ymin,ymax]); %%% %%% figure(fig_num); fig_num=fig_num+1; subplot(3,1,1); plot(freq,H2_p); out1=sprintf('H2 %s Frequency Response Function',tlabel); title(out1); grid on; ylabel('Phase (deg)'); axis([f1,f2,-180,180]); set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','lin','ytick',[-180,-90,0,90,180]); % subplot(3,1,[2 3]); plot(freq,H2_m); grid on; xlabel('Frequency(Hz)'); ylabel('G_X_X/G_X_F'); set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log',... 'YScale','log'); % ymax=10^(ceil(log10(max(H2_m)))); ymin=ymax; for i=1:mH if(H2_m(i)<ymin && freq(i)>f1 && freq(i)<f2 ) ymin=H2_m(i); end end % ymin=10^(floor(log10(ymin))); % if(ymin<ymax/10000) ymin=ymax/10000; end % axis([f1,f2,ymin,ymax]); %%% figure(fig_num); fig_num=fig_num+1; plot(tt,t_real); out1=sprintf('%s Impulse Response Function',tlabel); title(out1); xlabel('Time(sec)'); grid on; %%% figure(fig_num); fig_num=fig_num+1; plot(freq,COH); grid on; set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log',... 'YScale','lin'); xlabel('Frequency(Hz)'); ylabel('(\gamma_x_y)^2'); title('Coherence'); ymin=0.; ymax=1.2; axis([f1,f2,ymin,ymax]); %%% freq=fix_size(freq); H1_m=fix_size(H1_m); H1_p=fix_size(H1_p); % [xmax,fmax]=find_max([freq H1_m]); % out5 = sprintf('\n Peak occurs at %10.5g Hz ',fmax); disp(out5) % clear H1_frf_mp; clear H2_frf_mp; clear H1_frf_complex; clear H2_frf_complex; % disp(' '); disp(' Output arrays: '); disp(' magnitude & phase: H1_frf_mp & H2_frf_mp'); disp(' complex: H1_frf_complex & H2_frf_complex'); % H1_frf_mp=[freq H1_m H1_p]; H2_frf_mp=[freq H2_m H2_p]; H1_frf_complex=[freq real(H1C) imag(H1C)]; H2_frf_complex=[freq real(H2C) imag(H2C)]; % end