满足怎么样的条件，才能使得czt实现频率细化呢？

liuhuanjun

满足怎么样的条件，才能使得czt实现频率细化呢？
是不是变换要在单位园上？
还要满足什么条件吗
很迷惑阿
请哪位提示以下

[ 本帖最后由 yejet 于 2006-8-31 20:26 编辑 ]

yangzj

参数中w=exp(-j*2*pi/M)   中M大于序列x的长度N就是频率细化了
因为直接做fft的频率分辨率是否fs/N,
CZT的频率分辨率就是间隔w,也就是fs/M.
不过这种细化不能消除已经出现的干涉现象,要想减小干涉就要用ZOOM FFT了

liuhuanjun

请问，czt有什么缺点吗？
在什么情况下，不可以使用呢

liuhuanjun

是不是czt本身能够产生干涉阿

yangzj

不是本身产生干涉，而是不能消除信号本身靠的很近的频率成分的干涉，也就是说它无法分辨出两个靠的很近的频率，因为它并不是真正意义上的提高了频率分辨率．

songzy41

我完全同意yangzj 的观点，CZT不能使频谱细化，要频谱细化还得用ZoomFFT。我在资源下载中心内上传了一篇《两种快速频域细化分析方法的研究》，其中作者对CZT和ZoomFFT作了详细的讨论、仿真计算和两者比较，可看到在频率细化上ZoomFFT优于CZT。

w89986581

1. function zfft(x,fs,fi,fa,nfft)
   
2. %x信号;sf采样率;fi下限频率;fa上限频率
   
3. %%%%%%%%%%%%%%%%%细化fft%%%%%%%%%%%%%%%%
   
4. 
   
5. sf=250000;
   
6. f1=48000;
   
7. f2=49000;
   
8. fi=45000;
   
9. fa=52000;
   
10. nfft=512*2;
    
11. t=0:1/sf:2;
    
12. x=sin(2*pi*f1*t)+2*sin(2*pi*f2*t)+1*randn(size(t));
    
13. 
    
14. np=floor(sf/(fa-fi)/2);
    
15. nt=length(x);                           %nt数据长度
    
16. nf=2^nextpow2(nt);
    
17. na=round(0.5*nf/np+1);
    
18. n=0:nt-1;
    
19. b=n*pi*(fi+fa)/sf;
    
20. y=x.*exp(-i*b);
    
21. b=fft(y,nf);
    
22. a(1:na)=b(1:na);
    
23. a(nf-na+1:nf)=b(nf-na+1:nf);
    
24. b=ifft(a,nf);                           %
    
25. c=b(1:np:nt);                           %
    
26. a=fft(c,nfft)*2/nfft;                   %
    
27. y2=zeros(1,nfft/2);                     %y2为细化功率谱
    
28. y2(1:nfft/4)=a(nfft-nfft/4+1:nfft);
    
29. y2(nfft/4+1:nfft/2)=a(1:nfft/4);        
    
30. n=0:(nfft/2-1);
    
31. f2=fi+n*2*(fa-fi)/nfft;                 %f2为细化频率序列;fa为细化下限频率;fi为细化上限频率;
    
32. y1=fft(x,nfft)*2/nfft;                  %sf为采样频率;
    
33. f1=n*sf/nfft;
    
34. ni=round(fi*nfft/sf+1);
    
35. na=round(fa*nfft/sf+1);
    
36. %绘制输入时程曲线
    
37. subplot(2,1,1);
    
38. t=0:1/sf:(nt-1)/sf;
    
39. nn=1:length(t);
    
40. plot(t(nn),x(nn));
    
41. xlabel('时间(s)');
    
42. ylabel('幅值');
    
43. grid on;
    
44. %在相同的频率范围下绘制幅频曲线图
    
45. % subplot(2,1,2);
    
46. % nn=ni:na;
    
47. % plot(f1(nn),abs(y1(nn)),':',f2,abs(y2));
    
48. % xlabel('frequency(Hz)');
    
49. % ylabel('幅值');
    
50. % legend('普通','细化');
    
51. % Pxx=10*log10(abs(y2).^2*nfft/fs)
    
52. Pxx=10*log10(abs(y2).^2)
    
53. plot(f2,Pxx,'r')
    
54. grid on
    
55. hold on
    
56. psd(x,nfft,sf)
    
57. % %打开文件输出细化后的数据
    
58. % fid=fopen(fno,'w');
    
59. % for k=1:nfft/2
    
60. %   fprintf(fid,'%f%f%f\n',f2(k),real(y2(k)),imag(y2(k)));
    
61. % end
    
62. %
    
63. % status=fclose(fid);
    

复制代码

[ 本帖最后由 yejet 于 2006-8-31 20:27 编辑 ]

songzy41

也试着把《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;

miao7mijao

楼上的您好：
我把你写的程序运行了，
nt=length(x);     %计算读入数据长度
说x未定义，那我应该怎么办呢？

yangzj

原帖由 miao7mijao 于 2006-9-18 13:47 发表
楼上的您好：
我把你写的程序运行了，
nt=length(x);     %计算读入数据长度
说x未定义，那我应该怎么办呢？

不要直接运行他,做为一个函数调用

songzy41

你运行是应运行调试程序，而不是funcion。

miao7mijao

songzy41：你好！
看了你写的程序，其中几个地方不是很明白，我想问一下！谢谢先！
1） fa=fi+0.5 * fs/np;     %最大细化截止频率，这个不明白
2）na=round(0.5 * nf/np+1);     %确定细化带宽的数据长度
还有调试程序中的
x=y(2,:)也没弄明白；
就先问这几个，其实还有好多地方不是很懂呢，等待你的答复！

songzy41

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列数据。

miao7mijao

songzy41
谢谢你！以后还要麻烦你！

miao7mijao

w89986581：
这个公式
np=floor(sf/(fa-fi)/2);
是不是写错了，应该为np=floor(sf/2×(fa-fi))吧？