如何进行FFT频谱分析？

一、概述

FFT（Fast Fourier Transform）是一种用于对信号进行频域分析的常见方法。它的原理是将一个时间域的信号转换为一个频域的信号，从而将信号在不同频率上的表现展现出来。在实际应用中，FFT被广泛应用于音频、图像、通讯等领域。

Matlab中提供了许多用于FFT频谱分析的函数。其中最常用的是fft函数。下面将详细介绍如何在Matlab中实现FFT频谱分析。

二、FFT基本原理

FFT的基本原理是通过傅里叶变换将时域上的信号，转换到频域上的频谱图像。

在数字信号中，我们可以看做连续信号的抽样和量化后得到的。因此，可以用离散时间傅里叶变换（DTFT）将离散信号从时域转换到频域。

$$
X(e^{jomega})= sum_{n=-infty}^{infty}x[n]e^{-jomega n}
$$

其中，$x[n]$为离散时间的信号，$X(e^{jomega})$为频域信号。

本质上，DTFT是傅里叶变换的连续傅里叶变换形式，它是傅里叶变换在离散时间域上的表示。但由于采样率为有限，因此DTFT不能在计算机中实现。

因此，我们需要将DTFT变为离散傅里叶变换（DFT），它同样可以将时域上的信号，转换到频域上的频谱图像。

$$
X[k]= sum_{n=0}^{N-1}x[n]e^{-jfrac{2pi}{N}nk}, k=0,1,2, dots N-1
$$

其中，$x[n]$为时域离散信号，$X[k]$为离散频域信号。我们将$x[n]$看做离散信号的离散值，将$omega=frac{2pi}{N}k$看做是$omega$上的$kDeltaomega$，其中$Deltaomega=frac{2pi}{N}$，则离散傅里叶变换表示为

$$
X[k]= sum_{n=0}^{N-1}x[n]e^{-j n kDeltaomega}, k=0,1,2, dots N-1
$$

计算复杂度为$O(N^2)$。实际上，我们用FFT（快速傅里叶变换）来计算DFT。FFT的复杂度为$O(Nlog_2N)$，比DFT要快得多。

因此，fft函数实现FFT频谱分析的算法是：

1.将数据进行分组。

2.进行Butterfly操作，将数据进行排序，并计算第一组结果。

3.重复1和2过程，最后得到FFT频谱。

为了方便，Matlab中提供了fft函数来实现FFT频谱分析。

三、使用FFT频谱分析数据

通过fft函数，我们可以方便地进行FFT频谱分析。下面将介绍Matlab中用于FFT频谱分析的函数。

fft函数

fft函数的基本格式如下：

“`
Y = fft(X)
Y = fft(X,n)
Y = fft(X,[],dim)
“`

其中，X为要进行FFT频谱分析的序列，n为FFT变换长度（默认为X序列的长度），dim为进行FFT变换的维度（默认为第一个非单一的维度）。

Y为FFT变换后的结果。当X为向量时，Y将是X长度的复数向量。当X为矩阵时，Y将是沿着指定维度的FFT频域数据。

例如：

“`matlab
%生成一个长度为8的向量
>> x = linspace(0,2*pi,8);
>> y = sin(x);

%进行FFT变换
>> Y = fft(y);

%输出结果
>> abs(Y)
ans =

1 1 0.7071 0.7071 0.7071 0.7071 1.0000 1.0000
“`

上述示例中，我们生成了一个长度为8的正弦波，然后用fft函数对其进行FFT变换。最后输出了FFT变换后的结果。

fftshift函数

当我们进行FFT变换后，频率特性的直流分量（即直流成分）和高频成分位于FFT变换结果的两端，而重要的频率成分则位于中间。为了更好地观察FFT变换结果，我们可以使用fftshift函数将频率变换结果进行转换。

fftshift函数的基本形式如下：

“`
Y = fftshift(X)
Y = fftshift(X,dim)
“`

其中，X为要进行fftshift的序列，dim为进行fftshift的维度（默认为第一个非单一的维度）。

Y为fftshift变换后的结果。

例如：

“`matlab
%生成一个4×4的矩阵
>> X = magic(4)

X =

16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1

%对矩阵进行FFT变换
>> Y = fft2(X);

%对变换结果进行fftshift变换
>> Yshift = fftshift(Y);

%输出结果
>> abs(Y)
ans =

34.0000 6.4645 6.0000 3.5355
6.4645 6.0000 7.4645 6.0000
6.0000 7.4645 6.0000 6.4645
3.5355 6.0000 6.4645 34.0000

>> abs(Yshift)
ans =

6.0000 7.4645 6.0000 6.4645
3.5355 6.0000 6.4645 34.0000
6.4645 6.0000 7.4645 6.0000
34.0000 6.4645 6.0000 3.5355
“`

上述示例中，我们生成了一个4×4的矩阵，然后用fft2函数对其进行FFT变换，最后用fftshift函数对FFT变换结果进行变换。可以看出，fftshift变换将原本在FFT变换结果两端的直流分量和高频成分移动到了中间，更方便观察FFT变换结果。

四、Matlab示例

我们可以通过以下示例来更加深入地了解FFT频谱分析在Matlab中的应用。

示例1：生成正弦波并进行FFT变换

“`matlab
%生成正弦波
t = linspace(0,1,1000);
x = sin(2*pi*50*t)+0.5*sin(2*pi*150*t);

%进行FFT变换
N = length(x);
Y = fft(x);
P = abs(Y/N).^2;
f = linspace(0,1,N)*1000;

%绘制频谱图
f1 = figure;
set(f1,’Color’,[1 1 1]);
plot(f,P);
xlabel(‘Frequency (Hz)’);
ylabel(‘Power’);
“`

上述示例中，我们首先生成了一个包含50Hz和150Hz正弦波的信号。然后进行FFT变换，得到频域的频谱图像。最后用plot函数绘制出频谱图像。

示例2：使用窗函数平滑干扰

“`matlab
%生成正弦波
t = linspace(0,1,1000);
x = sin(2*pi*50*t)+0.5*sin(2*pi*150*t);

%加入干扰
xd = x + 2.5*randn(size(t));

%生成窗函数
win = hann(50);

%对干扰信号进行平滑
for i=1:950
xw(i) = mean(xd(i:i+49).*win’);
end

%进行FFT变换
N = length(xw);
Y = fft(xw);
P = abs(Y/N).^2;
f = linspace(0,1,N)*1000;

%绘制频谱图
f2 = figure;
set(f2,’Color’,[1 1 1]);
plot(f,P);
xlabel(‘Frequency (Hz)’);
ylabel(‘Power’);
“`

上述示例中，我们首先生成了一个包含50Hz和150Hz正弦波的信号，并加入高斯白噪声干扰。然后，通过hann窗函数进行平滑操作，滤除干扰信号，最后用FFT变换得到频域的频谱图像。可以看出，使用窗函数平滑干扰信号后，频谱图像清晰地显示了信号包含的频率成分。

示例3：实时进行FFT变换

“`matlab
%打开音频设备
recorder = audiorecorder(44100,16,1);
record(recorder);

%实时进行FFT变换
while(1)
if length(recorder) >= 44100
data = getaudiodata(recorder);
N = length(data);
Y = fft(data);
P = abs(Y/N).^2;
f = linspace(0,1,N)*44100;
plot(f,P);
xlabel(‘Frequency (Hz)’);
ylabel(‘Power’);
end
end
“`

上述示例中，我们使用Matlab内置的audiorecorder函数打开音频设备，然后在循环中使用getaudiodata函数获取实时音频信息。每当获取到1秒钟的音频数据后，就进行FFT变换并实时绘制频域的频谱图像。

通过上述示例，我们可以看到FFT频谱分析在音频处理中的应用。更多应用，等着程序员们去延伸。