基于matlab的信号频域的相干分析

% Copyright by Xiaoliang Meng % % Coded on August 1st, 2011 % % Modified on August 16th, 2011 % % Calculate the spectrum of the amplitude and phase angle of the cross power% % spectrum of two signals % % values returned: % % f ------ frequency % % z ------ cross power spectrum density % function myphase clear; wn = pwd; dn = wn; [filename,pathname] = uigetfile(dn,'请选择数据'); cd(wn); prompt = {'Total number of sensors:','First channel to do the CrossSpectral analysis:','Second channel to do the CrossSpectral analysis:'}; def = {'5','1','2'}; dlgTitle = 'Input the channel numbers'; lineNo = 1; answer = inputdlg(prompt,dlgTitle,lineNo,def); ChannelNum =str2num(answer{1}); SigNum1 =str2num(answer{2}); SigNum2 =str2num(answer{3}); DataCloumnNum = ChannelNum+1; SigNum1 = SigNum1+1; SigNum2 = SigNum2+1; fn = strcat(pathname,filename); fid = fopen(fn,'r'); % for k = 1:6 % tline = fgetl(fid); % end tline = fgetl(fid); % DISP = fscanf(fid,'%g %g %g %g %g %g %g',[7,inf]); % DISP = fscanf(fid,'%g %g %g %g',[4,inf]); DISP = fscanf(fid,'%g',[DataCloumnNum,inf]); fclose(fid); DISP = DISP'; Time = DISP(:,1); [m,n] = size(DISP); %去均值 Dm = mean(DISP(:,2:n)); for k = 2:n DISP(:,k) = DISP(:,k)-Dm(k-1); end sf = 1.0/(Time(2)-Time(1)); nfft = 2^(nextpow2(m)-1); % nfft = 4096; %建立频率向量 f = 0:sf/nfft:sf/2; f = f'; w = hanning(nfft); z = csd(DISP(:,SigNum1),DISP(:,SigNum2),nfft,sf,w,nfft/2); %绘制两个互谱信号的时程曲线图 figure(1); subplot(2,1,1); plot(Time,DISP(:,SigNum1)); xlabel('时间(s)'); ylabel('幅值'); grid on; title('信号时程曲线','Units','normalized ','Position',[0.5 1],'FontSize',13); subplot(2,1,2); plot(Time,DISP(:,SigNum2)); xlabel('时间(s)'); ylabel('幅值'); grid on; %绘制幅频曲线图 [maxvalue,maxindex] = max(abs(z)); maxfreq = f(maxindex); phase = angle(z(maxindex)); phasedeg = rad2deg(phase); figure(2); nn = 1:nfft/8; subplot(2,1,1); plot(f(nn),abs(z(nn)),'b-',maxfreq,maxvalue,'ro'); xlabel('频率(Hz)'); ylabel('幅值'); textlable=strcat('\leftarrow f= ',num2str(maxfreq)); text(maxfreq,maxvalue,textlable,'Fontsize',14,'Color','r'); grid on; title('互谱幅频、相频曲线','Units','normalized ','Position',[0.5 1],'FontSize',13); subplot(2,1,2); plot(f(nn),angle(z(nn)),'b-',maxfreq,phase,'ro'); xlabel('频率(Hz)'); ylabel('相位'); textlable=strcat('\leftarrow phase= ',num2str(phase)); text(maxfreq,phase,textlable,'Fontsize',14,'Color','r'); grid on; disp(['**** The phase between ',answer{2},' and ',answer{3},' is: ',num2str(phasedeg),' deg. ****']); timedelay = phase/(2*pi*maxfreq); disp(['**** The delaytime between ',answer{2},' and ',answer{3},' is: ',num2str(timedelay*1000),' ms. ****']); outputfn = strcat(pathname,'CrossSpectrum.dat'); fid = fopen(outputfn,'w'); fprintf(fid,' 频率(Hz) \t实部 \t虚部 \t幅值 \t相位(rad)\n'); for k =1:nfft/2 fprintf(fid,'%12.6f\t%12.4e\t%12.4e\t%12.4e\t%12.4e\n',f(k),real(z(k)),imag(z(k)),abs(z(k)),angle(z(k))); end fclose(fid);