[原创]exzfft_m

songzy41

在http://forum.vibunion.com/forum/vi ... 23569&highlight=czt上我曾提供了zfft_m的函数，yangzj 指出：“这个程序的采用这样的滤波有问题哦,使整个程序都没有意义了.既然都已经做了 nf点的FFT了,那直接用这个结果的频率分辨率已经很低了,也就没必要后面的细化了. ”他批评得对。我把原函数作了改进，用了decimate函数既滤波又下采样，现提供该ZFFT函数共享。该函数文件名为exzfft_m:
function y=exzfft_m(x,fi,fs,nfft,np)
% x 被测信号，被测信号长度要>=nfft*np
% 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);
%  滤波和下采样
c=decimate(y,np);
% FFT
y=fft(c, nfft) * 2/nfft;
以下又给出测试程序，计算结果如下：

[ 本帖最后由 zhangnan3509 于 2007-7-4 15:32 编辑 ]

hnlzx

运行正常，辛苦了，谢谢！

xiaodinger

测试程序怎么不见了？还能提供吗？:loveliness:

songzy41

附上测试程序
%
% test_exzfftm
clear all; clc;

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=exzfft_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;
并附上数据文件

[ 本帖最后由 songzy41 于 2007-9-20 09:35 编辑 ]

xiaodinger

就上述程序而言，如果要求出细化后的频率和所在的谱线位置，该怎么做呢？我的目的是看做完zfft后，精度能达到多少？

xiaodinger

程序运行时出现提示：
??? Input argument 'x' is undefined.

Error in ==> C:\MATLAB6p5p1\work\exzfft_m.m
On line 8  ==> nt=length(x);
是不是x=y(2,:);有问题呀？

songzy41

1，我又试了，寿行没有问题；
2，“要求出细化后的频率和所在的谱线位置”，基本上和FFT计算后的求法一样。但由于频移了，所以y2(1)对应的频率不是0频率，而是最小细化截止频率fi，而频率的计算就是f2的式子。
这给出的仅是频谱的细化程序，即提高了频谱的分辨率；如果楼主想精确地求出最大值的幅值和频率，还是使用修正法较好。

xiaodinger

请问fi和fa确定的依据是什么？如果信号的频率未知（但是有一个范围），那么该怎样确定细化的最低频率和最大截至频率呢？难道实际运行时，fi都必须是预先给定吗？如果最大截止频率fa=fi+0.5 * sf/np，那么改变np时，fa就可能会比所测频率f小了。

songzy41

可以参看一下帖子http://forum.vibunion.com/forum/vi ... 23569&highlight=czt，我这里的exzfft_m是从该贴子上的zfft_m修改而来的。而zfft_m又是从王济编写的《MATLAB在振动信号处理中的应用》的一程序改来的。在王济编写的程序中，fi和fa都是设定的，所以不存在fa可能会比所测频率f小的问题。
同样在王济的程序中，它是把整个数据都做了高分辨率的FFT变换，正好适合楼主不知信号的频率的分布的情形，楼主可上这个帖子看一下。

witty01

四楼里面
y2(1: nfft/4) =a(nfft-nfft/4+1 : nfft);
y2(nfft/4+1 : nfft/2) =a(1: nfft/4);

这两句有什么依据呢，没想通

songzy41

原帖由 witty01 于 2007-10-31 12:48 发表
四楼里面
y2(1: nfft/4) =a(nfft-nfft/4+1 : nfft);
y2(nfft/4+1 : nfft/2) =a(1: nfft/4);

这两句有什么依据呢，没想通

在exzfft_m中我们可以看到是进行了频移，把(fi+fa)/2移到了0频率，所以返回后把负频率部分移回到正频部分，把谱线重新安排，以上的两语句就是做这一工作的。

witty01

感谢songzy41的回答，我想我还得再问一下，

在把谱线重新排列时，把zfft的结果平均分成了四部分，为什么取了第四部分和第一部分分别付给最后的y2的前半段和后半段呢？

fft的结果是对称的，显示时只要前面一半就可以了，这个为什么不是？谢谢

yangzj

移频后是一个复信号，不再是对称的了。你把每一步的频域变化过程用一个图来表达，就很清楚了。

songzy41

原帖由 witty01 于 2007-11-5 16:30 发表
感谢songzy41的回答，我想我还得再问一下，

在把谱线重新排列时，把zfft的结果平均分成了四部分，为什么取了第四部分和第一部分分别付给最后的y2的前半段和后半段呢？

fft的结果是对称的，显示时只要前面 ...

大版主yangzj 说得对。FFT对称是有条件的，即实数序列的FFT将是实部偶对称，虚部奇对称。但在exzfft_m中，为了细化，把输入序列x乘以旋转因子进行频移：
b=n*pi* (fi+fa)/fs;     %乘单位旋转因子进行频移
y=x.*exp(-i*b);
y便是一个复数序列了，即使经过滤波和下采样，依然是一个复数序列，y的FFT便是复数序列的FFT，FFT后不再满足对称的性质了。
同样在细化以后，采样频率为sf/np，信号复调 制后0频率为原信号的 (fi+fa)/2，信号一样满足采样定理，信号的最高频率应小于sf/np/2，对应于原信号便是在fi~fa之间，即对应于y2的第四部分和第一部分。

witty01

:lol ,解释的很细致，谢谢你们！！