Matlab实现FFT变换

文章目录

• • Matlab实现FFT变换
  • • 原理
    • 实现
    • • 手算验证
      • 简单fft变换和频谱
      • 求取功率谱
    • 结论

原理

FFT是一种计算离散傅里叶变换（DFT）的快速算法。DFT将时域上的信号转换为频域上的信号，可以用以下公式表示：

Xk=∑n=0N−1xne−i2πkn/NX_k=\sum_{n=0}^{N-1}x_n e^{-i2\pi kn/N}Xk​=n=0∑N−1​xn​e−i2πkn/N

其中，xnx_nxn​是时域上的信号序列，XkX_kXk​是频域上的信号序列，kkk为频率编号（0≤k

FFT算法通过分治策略将DFT算法的计算复杂度从O(N2)O(N^2)O(N2)降低到O(Nlog2N)O(Nlog_2N)O(Nlog2​N)，从而实现了在计算机上快速计算DFT的目的。

实现

手算验证

在这里，我们将给出一个简单的例子来说明如何使用Matlab进行FFT变换。我们首先生成一个简单的数组：

x = [1, 2, 3, 4, 5, 6, 7, 8];

接下来，我们使用Matlab内置的fft函数对这个数组进行FFT变换：

X = fft(x);

这个操作会返回一个和输入数组长度相同的复数数组。我们可以使用Matlab的disp函数打印出这个数组：

disp(X);

这个例子的输出结果如下：

   36.0000 + 0.0000i-4.0000 + 9.6569i-4.0000 + 4.0000i-4.0000 + 1.6569i-4.0000 + 0.0000i-4.0000 - 1.6569i-4.0000 - 4.0000i-4.0000 - 9.6569i

可以看到，输出结果是一个长度为8的复数数组。

为了验证FFT的正确性，我们可以手动计算这个输入数组的FFT结果，然后将结果与Matlab计算的结果进行比较。FFT算法的计算过程可以用以下公式表示：

Xk=∑n=0N−1xne−i2πkn/NX_k=\sum_{n=0}^{N-1}x_n e^{-i2\pi kn/N}Xk​=n=0∑N−1​xn​e−i2πkn/N

其中，xnx_nxn​是时域上的信号序列，XkX_kXk​是频域上的信号序列，kkk为频率编号（0≤k

对于输入数组x=[1,2,3,4,5,6,7,8]x=[1,2,3,4,5,6,7,8]x=[1,2,3,4,5,6,7,8]，我们有N=8N=8N=8。因此，X0X_0X0​的计算公式为：

X0=∑n=07xne−i2π0n/8=36X_0 = \sum_{n=0}^{7}x_n e^{-i2\pi 0n/8}=36X0​=n=0∑7​xn​e−i2π0n/8=36

接下来，我们可以计算X1X_1X1​：

X1=∑n=07xne−i2π1n/8=−4+9.6569iX_1 = \sum_{n=0}^{7}x_n e^{-i2\pi 1n/8}=-4+9.6569iX1​=n=0∑7​xn​e−i2π1n/8=−4+9.6569i

以此类推，我们可以计算出所有的XkX_kXk​。最终结果应该与Matlab计算的结果一致。

简单fft变换和频谱

以下是一个简单的Matlab代码实现FFT变换：

% 生成测试信号
Fs = 1000;      % 采样频率
t = 0:1/Fs:1-1/Fs;   % 时间向量
x = 1*sin(2*pi*100*t); % 信号% 绘制信号图
subplot(2,1,1);
plot(t,x);
title('信号');
xlabel('时间 (s)');
ylabel('幅度');% 计算FFT
N = length(x);
X = fft(x);
f = Fs*(0:(N/2))/N;% 绘制FFT图
subplot(2,1,2);
plot(f,abs(X(1:N/2+1)));
title('FFT');
xlabel('频率 (Hz)');
ylabel('幅度');

在这个例子中，我们生成了一个频率为100Hz的正弦信号。我们使用Matlab的fft函数计算FFT，并将结果绘制成幅度谱。注意，在绘制幅度谱时，我们只绘制了频率为正的一半，因为FFT算法输出的结果是对称的。

求取功率谱

通过FFT变换可以得到信号的幅度谱，但是为了更好地了解信号特性，我们通常需要求取信号的功率谱密度。功率谱密度描述了信号在不同频率下的功率分布情况。

求取功率谱的方法是，将信号进行FFT变换后，将每个频率上的幅度平方除以信号长度，并乘以一个系数，即可得到功率谱密度。具体公式如下：

Pk=2∣Xk∣2NP_k=\frac{2|X_k|^2}{N}Pk​=N2∣Xk​∣2​

其中，PkP_kPk​是频率为kkk的功率谱密度，XkX_kXk​是频率为kkk的信号幅度，NNN是信号长度。

以下是一个简单的Matlab代码实现求取功率谱：

% 生成测试信号
Fs = 1000;      % 采样频率
t = 0:1/Fs:1-1/Fs;   % 时间向量
x = 1*sin(2*pi*100*t); % 信号% 绘制信号图
subplot(2,1,1);
plot(t,x);
title('信号');
xlabel('时间 (s)');
ylabel('幅度');% 计算FFT
N = length(x);
X = fft(x);
f = Fs*(0:(N/2))/N;% 计算功率谱
P = (2*abs(X(1:N/2+1)).^2)/N;% 绘制功率谱图
subplot(2,1,2);
plot(f,P);
title('功率谱密度');
xlabel('频率 (Hz)');
ylabel('功率');

注意，我们在计算功率谱时，使用了一个系数222，这是因为我们只绘制了频率为正的一半，而实际上信号的功率谱是对称的。

结论

本文介绍了如何使用Matlab实现FFT变换，并求取信号的功率谱密度。通过FFT变换，我们可以将信号从时域转换到频域，进一步了解信号的特性。