频谱细化-----FFT+FT算法介绍及MATLAB实现

FFT+FT算法原理

设信号的长度为$M$，设FFT点数为$N$，$N>M$，即对信号补零进行FFT，$N$点FFT算法对应的DFT变换公式为
$$X\left(k\right)=\sum _{n=0}^{N-1}x\left(n\right)W_N^{kn},k=0,1,...,N-1$$
经过补零后的$N$点FFT算法对应的DFT变换公式为
$$X\left(k\right)=\sum _{n=0}^{N-1}x\left(n\right)e^{-j\frac{2\pi }{N}kn}$$
由于补零后的信号从$M$到$N-1$之间都是0，那么
$$X\left(k\right)=\sum _{n=0}^{M-1}x\left(n\right)e^{-j\frac{2\pi }{N}kn}$$
设从FFT算法计算出的频谱中寻找到频谱最大值对应谱线的标号为$k_0$,次大值频谱对应的谱线标号为$k_1$，即
$$X\left(k_0\right)=\sum _{n=0}^{M-1}x\left(n\right)e^{-j2\pi \frac{k_0}{N}n}$$
$$X\left(k_1\right)=\sum _{n=0}^{M-1}x\left(n\right)e^{-j2\pi \frac{k_1}{N}n}$$
由于栅栏效应导致频谱最大值与频谱次大值之间的频率没有被检测出来，为了提高精度需要得到$k_0/N$和$k_1/M=\left(k_0\pm 1\right)/M$更多的值，也就是对初期通过FFT算法计算得到的频谱做进一步的细化处理。
设在原相邻频谱线之间再细化$K$倍，也就是需要求得频谱在如下频点的频谱幅值：
$$\frac{K\cdot k_0}{K\cdot N},\frac{K\cdot k_0+1}{K\cdot N},...,\frac{K\cdot k_0+K-1}{K\cdot N},\frac{K\cdot k_0+K}{K\cdot N},k_1=k_0+1$$
或在如下频点的频谱幅值
$$\frac{K\cdot k_1}{K\cdot N},\frac{K\cdot k_1+1}{K\cdot N},...,\frac{K\cdot k_1+K-1}{K\cdot N},\frac{K\cdot k_1+K}{K\cdot N},k_1=k_0-1$$
计算上述频点对应的频谱的时候，只需要调用DFT公式计算，不需要使用DFT计算整个频谱。
下面对上述算法进行仿真
参数设置为

参数名称	参数值
信号长度	2048
采样率	4096
FFT点数	2048
信号频率	346.6Hz

仿真结果如下：

从图中可以看出，直接进行FFT变换的话是无法分辨出频率为346.6Hz的信号，但是经过FFT+FT算法细化之后可以得到信号的信号的频率的，也体现了频谱细化的好处。
代码如下：

clear;
close all;
clc;
N=1024;                 %信号长度
fs=4096;                %采样率
t=(0:N-1)/fs;  
x=exp(-2*pi*0.005*346.6.*t).*cos(2*pi*346.6.*t+pi/3); %信号
y=fft(x,1024); %信号傅氏变换
figure(1) 
f1=(0:N-1)/N*fs; 
figure(1);
subplot(211);
plot(f1(1:N/2),2/N*abs(y(1:N/2)));
xlabel('f/Hz');
title('直接进行FFT')
f=342:0.1:352;    
for k=1:length(f)
    xf(k)=2/N*sum(x.*exp(-1j*2*pi*(0:N-1)*f(k)/fs)); 
end     
subplot(212);
plot(f,abs(xf)) 
title('FFT+FT算法处理后')