|
ZFFT程序(Frome《MATLAB在振动信号处理中的应用》)
- function zfft(x,fs,fi,fa,nfft)
- %x信号;sf采样率;fi下限频率;fa上限频率
- %%%%%%%%%%%%%%%%%细化fft%%%%%%%%%%%%%%%%
- sf=250000;
- f1=48000;
- f2=49000;
- fi=45000;
- fa=52000;
- nfft=512*2;
- t=0:1/sf:2;
- x=sin(2*pi*f1*t)+2*sin(2*pi*f2*t)+1*randn(size(t));
- np=floor(sf/(fa-fi)/2);
- nt=length(x); %nt数据长度
- nf=2^nextpow2(nt);
- na=round(0.5*nf/np+1);
- n=0:nt-1;
- b=n*pi*(fi+fa)/sf;
- y=x.*exp(-i*b);
- b=fft(y,nf);
- a(1:na)=b(1:na);
- a(nf-na+1:nf)=b(nf-na+1:nf);
- b=ifft(a,nf); %
- c=b(1:np:nt); %
- a=fft(c,nfft)*2/nfft; %
- y2=zeros(1,nfft/2); %y2为细化功率谱
- y2(1:nfft/4)=a(nfft-nfft/4+1:nfft);
- y2(nfft/4+1:nfft/2)=a(1:nfft/4);
- n=0:(nfft/2-1);
- f2=fi+n*2*(fa-fi)/nfft; %f2为细化频率序列;fa为细化下限频率;fi为细化上限频率;
- y1=fft(x,nfft)*2/nfft; %sf为采样频率;
- f1=n*sf/nfft;
- ni=round(fi*nfft/sf+1);
- na=round(fa*nfft/sf+1);
- %绘制输入时程曲线
- subplot(2,1,1);
- t=0:1/sf:(nt-1)/sf;
- nn=1:length(t);
- plot(t(nn),x(nn));
- xlabel('时间(s)');
- ylabel('幅值');
- grid on;
- %在相同的频率范围下绘制幅频曲线图
- % subplot(2,1,2);
- % nn=ni:na;
- % plot(f1(nn),abs(y1(nn)),':',f2,abs(y2));
- % xlabel('frequency(Hz)');
- % ylabel('幅值');
- % legend('普通','细化');
- % Pxx=10*log10(abs(y2).^2*nfft/fs)
- Pxx=10*log10(abs(y2).^2)
- plot(f2,Pxx,'r')
- grid on
- hold on
- psd(x,nfft,sf)
- % %打开文件输出细化后的数据
- % fid=fopen(fno,'w');
- % for k=1:nfft/2
- % fprintf(fid,'%f%f%f\n',f2(k),real(y2(k)),imag(y2(k)));
- % end
- %
- % status=fclose(fid);
复制代码
[ 本帖最后由 yejet 于 2006-8-31 20:27 编辑 ] |
评分
-
1
查看全部评分
-
|