声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 9638|回复: 28

[FFT] 满足怎么样的条件,才能使得czt实现频率细化呢?

[复制链接]
发表于 2006-8-28 09:59 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

[ 本帖最后由 yejet 于 2006-8-31 20:26 编辑 ]
回复
分享到:

使用道具 举报

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

评分

1

查看全部评分

 楼主| 发表于 2006-8-28 12:01 | 显示全部楼层
请问,czt有什么缺点吗?
在什么情况下,不可以使用呢
 楼主| 发表于 2006-8-28 12:05 | 显示全部楼层
是不是czt本身能够产生干涉阿
发表于 2006-8-28 12:47 | 显示全部楼层
不是本身产生干涉,而是不能消除信号本身靠的很近的频率成分的干涉,也就是说它无法分辨出两个靠的很近的频率,因为它并不是真正意义上的提高了频率分辨率.
发表于 2006-8-28 16:11 | 显示全部楼层
我完全同意yangzj 的观点,CZT不能使频谱细化,要频谱细化还得用ZoomFFT。我在资源下载中心内上传了一篇《两种快速频域细化分析方法的研究》,其中作者对CZT和ZoomFFT作了详细的讨论、仿真计算和两者比较,可看到在频率细化上ZoomFFT优于CZT。

评分

1

查看全部评分

发表于 2006-8-29 16:26 | 显示全部楼层

ZFFT程序(Frome《MATLAB在振动信号处理中的应用》)

  1. function zfft(x,fs,fi,fa,nfft)
  2. %x信号;sf采样率;fi下限频率;fa上限频率
  3. %%%%%%%%%%%%%%%%%细化fft%%%%%%%%%%%%%%%%

  4. sf=250000;
  5. f1=48000;
  6. f2=49000;
  7. fi=45000;
  8. fa=52000;
  9. nfft=512*2;
  10. t=0:1/sf:2;
  11. x=sin(2*pi*f1*t)+2*sin(2*pi*f2*t)+1*randn(size(t));

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

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

评分

1

查看全部评分

发表于 2006-8-31 16:44 | 显示全部楼层
也试着把《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;

图谱比较

图谱比较

xdata.txt

213.03 KB, 下载次数: 215

数据

评分

1

查看全部评分

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



不要直接运行他,做为一个函数调用
发表于 2006-9-18 16:06 | 显示全部楼层
你运行是应运行调试程序,而不是funcion。
发表于 2006-9-19 13:56 | 显示全部楼层
songzy41:你好!
看了你写的程序,其中几个地方不是很明白,我想问一下!谢谢先!
1) fa=fi+0.5 * fs/np;     %最大细化截止频率,这个不明白
2)na=round(0.5 * nf/np+1);     %确定细化带宽的数据长度
还有调试程序中的
x=y(2,:)也没弄明白;
就先问这几个,其实还有好多地方不是很懂呢,等待你的答复!
发表于 2006-9-19 15:33 | 显示全部楼层
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列数据。
发表于 2006-9-19 16:48 | 显示全部楼层
songzy41
谢谢你!以后还要麻烦你!
发表于 2006-9-20 09:33 | 显示全部楼层
w89986581:
这个公式
np=floor(sf/(fa-fi)/2);
是不是写错了,应该为np=floor(sf/2×(fa-fi))吧?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 20:45 , Processed in 0.187903 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表