满足怎么样的条件,才能使得czt实现频率细化呢?
满足怎么样的条件,才能使得czt实现频率细化呢?是不是变换要在单位园上?
还要满足什么条件吗
很迷惑阿
请哪位提示以下
[ 本帖最后由 yejet 于 2006-8-31 20:26 编辑 ] 参数中w=exp(-j*2*pi/M) 中M大于序列x的长度N就是频率细化了
因为直接做fft的频率分辨率是否fs/N,
CZT的频率分辨率就是间隔w,也就是fs/M.
不过这种细化不能消除已经出现的干涉现象,要想减小干涉就要用ZOOM FFT了 请问,czt有什么缺点吗?
在什么情况下,不可以使用呢 是不是czt本身能够产生干涉阿 不是本身产生干涉,而是不能消除信号本身靠的很近的频率成分的干涉,也就是说它无法分辨出两个靠的很近的频率,因为它并不是真正意义上的提高了频率分辨率. 我完全同意yangzj 的观点,CZT不能使频谱细化,要频谱细化还得用ZoomFFT。我在资源下载中心内上传了一篇《两种快速频域细化分析方法的研究》,其中作者对CZT和ZoomFFT作了详细的讨论、仿真计算和两者比较,可看到在频率细化上ZoomFFT优于CZT。
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 编辑 ] 也试着把《MATLAB在振动信号处理中的应用》中的ZFFT改为函数,方便其它坛友的调用,并提供调试程序和数据。下图是调试程序给出的计算结果。
函数为:
function y=zfft_m(x,fi,fs,nfft,np)
% x 被测信号
% fi细化的最低频率
% fs采样频率
% nfft 作细化FFT长
% np 放大倍数
% y细化FFT输出
nt=length(x); %计算读入数据长度
fa=fi+0.5 * fs/np; %最大细化截止频率
nf= 2^nextpow2(nt); %取大于nt且最接近nt的整数次方为FFT长度
na=round(0.5 * nf/np+1); %确定细化带宽的数据长度
% 频移
n=0: nt-1; %建一个递增向量
b=n*pi* (fi+fa)/fs; %乘单位旋转因子进行频移
y=x.*exp(-i*b);
%FFT
b= fft(y, nf);
% 低通滤波
a(1: na) =b(1: na);
a(nf-na+1 : nf) =b(nf-na+1 : nf);
% IFFT
b=ifft(a, nf);
%下采样
c=b(1 : np: nt);
% 再一次FFT
y=fft(c, nfft) * 2/nfft;
调试程序:
sf=200; %采样频率
fi=5; % 最小细化截止频率
np=8; % 放大倍数
nfft=1024; % FFT长度
y=load('xdata.txt'); %读入数据
y=y';
x=y(2,:);
fa=fi+0.5 * sf/np;
nt=length(x);
% 计算zfft
a=zfft_m(x,fi,sf,nfft,np);
% 排列数据
y2=zeros(1, nfft/2);
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;
% 把信号按nfft长作FFT计算
y1=fft(x, nfft) * 2/nfft;
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 : 3000;
plot (t(nn), x(nn));
xlabel ('Time (s)');
ylabel ('Amplitude');
title('Waveform');
grid on;
subplot (2, 1, 2);
nn= ni :na;
plot (f1(nn), abs(y1(nn)),':',f2, abs(y2));
xlabel ('Frequency (Hz)');
ylabel ('Amplitude');
legend ('未细化' ,'细化');
grid on; 楼上的您好:
我把你写的程序运行了,
nt=length(x); %计算读入数据长度
说x未定义,那我应该怎么办呢? 原帖由 miao7mijao 于 2006-9-18 13:47 发表
楼上的您好:
我把你写的程序运行了,
nt=length(x); %计算读入数据长度
说x未定义,那我应该怎么办呢?
不要直接运行他,做为一个函数调用 你运行是应运行调试程序,而不是funcion。 songzy41:你好!
看了你写的程序,其中几个地方不是很明白,我想问一下!谢谢先!
1) fa=fi+0.5 * fs/np; %最大细化截止频率,这个不明白
2)na=round(0.5 * nf/np+1); %确定细化带宽的数据长度
还有调试程序中的
x=y(2,:)也没弄明白;
就先问这几个,其实还有好多地方不是很懂呢,等待你的答复! 1,当采样频率为fs时,在频谱中最高可显示fs/2=0.5*fs。现细化放大是np倍,最高显示减小到0.5*fs/np;在细化时经频移,fi移到了0频,这样最大细化截止频率就为fa=fi+0.5 * fs/np。
2,首先要搞明白na是用在什么地方的。在程序中na是用于低通滤波,即在数据x经FFT后,只取0~na之间的谱线,而na以外的谱线(在正频率区间)都置0。
nf是x数据长nt最接近的2的幂次,以后的FFT便按nf来计算。计算出来有nf条谱线,正频率部分只需取其中的一半,即0.5 *nf条谱线。细化放大倍数为np,因此细化的区间仅在0~0.5 *nf/np之间。但0.5 *nf/np可能不是整数,谱线的具体条数必须是整数,故用了na=round(0.5 * nf/np+1);
3,读入后的y是有2列的数组,笫1列代表时间,笫2列代表测量到的振动数据。x=y(2,:)表示取笫2列数据。 songzy41
谢谢你!以后还要麻烦你! w89986581:
这个公式
np=floor(sf/(fa-fi)/2);
是不是写错了,应该为np=floor(sf/2×(fa-fi))吧?
页:
[1]
2