傅里叶变换

傅里叶变换

傅里叶变换是一个数学公式，用于将按时间或空间采样的信号变换为按时序或空间频率采样的相同信号。在信号处理中，傅里叶变换可以揭示信号的重要特征（即其频率分量）。

对于包含 n 个均匀采样点的向量 x，其傅里叶变换定义为

MATLAB® 中的 fft 函数使用快速傅里叶变换算法来计算数据的傅里叶变换。以正弦信号 x 为例，该信号是时间 t 的函数，频率分量为 15 Hz 和 20 Hz。使用在 10 秒周期内以 1/50 秒为增量进行采样的时间向量。

Ts = 1/50;
t = 0:Ts:10-Ts;
x = sin(2*pi*15*t) + sin(2*pi*20*t);
plot(t,x)
xlabel('Time (seconds)')
ylabel('Amplitude')

计算信号的傅里叶变换，并在频率空间创建对应于信号采样的向量 f。

y = fft(x);
fs = 1/Ts;
f = (0:length(y)-1)*fs/length(y);

以频率函数形式绘制信号幅值时，幅值尖峰对应于信号的 15 Hz 和 20 Hz 频率分量。

plot(f,abs(y))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Magnitude')

该变换还会生成尖峰的镜像副本，该副本对应于信号的负频率。为了更好地以可视化方式呈现周期性，您可以使用 fftshift 函数对变换执行以零为中心的循环平移。

n = length(x);
fshift = (-n/2:n/2-1)*(fs/n);
yshift = fftshift(y);
plot(fshift,abs(yshift))
xlabel('Frequency (Hz)')
ylabel('Magnitude')

含噪信号

在科学应用中，信号经常遭到随机噪声破坏，掩盖其频率分量。傅里叶变换可以清除随机噪声并显现频率。例如，通过在原始信号 x 中注入高斯噪声，创建一个新信号 xnoise。

rng('default')
xnoise = x + 2.5*randn(size(t));

频率函数形式的信号功率是信号处理中的一种常用度量。功率是信号的傅里叶变换按频率样本数进行归一化后的平方幅值。计算并绘制以零频率为中心的含噪信号的功率谱。尽管存在噪声，您仍可以根据功率中的尖峰辨识出信号的频率。

ynoise = fft(xnoise);
ynoiseshift = fftshift(ynoise);
power = abs(ynoiseshift).^2/n; 
plot(fshift,power)
title('Power')
xlabel('Frequency (Hz)')
ylabel('Power')

计算效率

直接使用傅里叶变换公式分别计算 y 的 n 个元素需要 n2 数量级的浮点运算。使用快速傅里叶变换算法，则只需要 n log n 数量级的运算。在处理包含成百上千万个数据点的数据时，这一计算效率会带来很大的优势。在 n 具有较小的质因数（例如 n 为 2 的幂）时，许多专门的快速傅里叶变换算法实现可进一步提高效率。

以加利福尼亚海岸的水下麦克风所收集的音频数据为例。在康奈尔大学生物声学研究项目维护的库中可以找到这些数据。载入包含太平洋蓝鲸鸣声的文件 bluewhale.au，并对其中一部分数据进行格式化。由于蓝鲸的叫声是低频声音，人类几乎听不到。数据中的时间标度压缩了 10 倍，以便提高音调并使叫声更清晰可闻。可使用命令 sound(x,fs) 来收听完整的音频文件。

whaleFile = 'bluewhale.au';
[x,fs] = audioread(whaleFile);
whaleMoan = x(2.45e4:3.10e4);
t = 10*(0:1/fs:(length(whaleMoan)-1)/fs);

plot(t,whaleMoan)
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 t(end)])

指定新的信号长度，该长度是大于原始长度的最邻近的 2 的幂。然后使用 fft 和新的信号长度计算傅里叶变换。fft 会自动用零填充数据，以增加样本大小。此填充操作可以大幅提高变换计算的速度，对于具有较大质因数的样本大小更是如此。

m = length(whaleMoan);
n = pow2(nextpow2(m));
y = fft(whaleMoan,n);

绘制信号的功率谱。绘图指示，呻吟音包含约 17 Hz 的基本频率和一系列谐波（其中强调了第二个谐波）。

f = (0:n-1)*(fs/n)/10; % frequency vector
power = abs(y).^2/n;   % power spectrum
plot(f(1:floor(n/2)),power(1:floor(n/2)))
xlabel('Frequency (Hz)')
ylabel('Power')

正弦波的相位

使用傅里叶变换，您还可以提取原始信号的相位频谱。例如，创建一个由频率为 15 Hz 和 40 Hz 的两个正弦波组成的信号。第一个正弦波是相位为 −π/4 的余弦波，第二个是相位为 π/2 的余弦波。以 100 Hz 的频率对信号采样 1 秒。

fs = 100;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*15*t - pi/4) - sin(2*pi*40*t);

计算信号的傅里叶变换。将变换幅值绘制为频率函数。

y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*fs;

stem(f,abs(z))
xlabel("Frequency (Hz)")
ylabel("|y|")
grid

计算变换的相位，删除小幅值变换值。将相位绘制为频率的函数。

tol = 1e-6;
z(abs(z) < tol) = 0;

theta = angle(z);

stem(f,theta/pi)
xlabel("Frequency (Hz)")
ylabel("Phase / \pi")
grid

离散傅里叶变换

         离散傅里叶变换（即 DFT）是数字信号处理的首要工具。该产品的基础是快速傅里叶变换 (FFT)，这是一种可减少执行时间的 DFT 计算方法。许多工具箱函数（包括 Z 域频率响应、频谱和倒频谱分析，以及一些滤波器设计和实现函数）都支持 FFT。

       MATLAB® 环境提供 fft 和 ifft 函数，分别用于计算离散傅里叶变换及其逆变换。对于输入序列 x 及其变换版本 X（围绕单位圆的等间距频率的离散时间傅里叶变换），这两个函数实现以下关系：

和

在这些方程中，序列下标从 1 而不是 0 开始，因为采用 MATLAB 向量索引方案，并且

注意：MATLAB 约定是对 fft 函数使用负 j。这是工程约定；物理和纯数学通常使用正 j。

       使用单个输入参数 x 的 fft 计算输入向量或矩阵的 DFT。如果 x 是向量，fft 计算向量的 DFT；如果 x 是矩形数组，fft 计算每个数组列的 DFT。

例如，创建时间向量和信号：

t = 0:1/100:10-1/100;                     % Time vector
x = sin(2*pi*15*t) + sin(2*pi*40*t);      % Signal

计算信号的 DFT 以及变换后的序列的幅值和相位。通过将小幅值变换值设置为零来减少计算相位时的舍入误差。

y = fft(x);                               % Compute DFT of x
m = abs(y);                               % Magnitude
y(m<1e-6) = 0;
p = unwrap(angle(y));                     % Phase

要以度为单位绘制幅值和相位，请键入以下命令：

f = (0:length(y)-1)*100/length(y);        % Frequency vector

subplot(2,1,1)
plot(f,m)
title('Magnitude')
ax = gca;
ax.XTick = [15 40 60 85];

subplot(2,1,2)
plot(f,p*180/pi)
title('Phase')
ax = gca;
ax.XTick = [15 40 60 85];

fft 的第二个参数指定变换的点数 n，表示 DFT 的长度：

n = 512;
y = fft(x,n);
m = abs(y);
p = unwrap(angle(y));
f = (0:length(y)-1)*100/length(y);

subplot(2,1,1)
plot(f,m)
title('Magnitude')
ax = gca;
ax.XTick = [15 40 60 85];

subplot(2,1,2)
plot(f,p*180/pi)
title('Phase')
ax = gca;
ax.XTick = [15 40 60 85];

       在本例中，如果输入序列比 n 短，fft 会用零填充输入序列，如果输入序列比 n 长，则会截断序列。如果未指定 n，则默认为输入序列的长度。fft 的执行时间取决于其执行的 DFT 的长度 n；有关该算法的详细信息，请参阅 fft 参考页。

注意：得到的 FFT 振幅是 A*n/2，其中 A 是原始振幅，n 是 FFT 点数。仅当 FFT 点的数量大于或等于数据样本的数量时，上述情形才成立。如果 FFT 点数小于数据样本数，则 FFT 振幅比原始振幅低上述量。

       离散傅里叶逆变换函数 ifft 也接受输入序列以及可选的变换所需点数。尝试以下示例；原始序列 x 和重新构造的序列是相同的（在舍入误差内）。

t = 0:1/255:1;
x = sin(2*pi*120*t);
y = real(ifft(fft(x)));

figure
plot(t,x-y)

       该工具箱还包括二维 FFT 及其逆变换的函数，即 fft2 和 ifft2。这些函数对于二维信号或图像处理非常有用。goertzel 函数是计算 DFT 的另一种算法，它也包含在工具箱中。此函数可高效计算长信号中一部分的 DFT。

        有时可以方便地重新排列 fft 或 fft2 函数的输出，使零频率分量位于序列的中心。函数 fftshift 将零频率分量移至向量或矩阵的中心。

另请参阅

fft | fftshift | nextpow2 | ifft | fft2 | fftn | fftw

相关主题

• 二维傅里叶变换

外部网站

• 傅里叶分析（MathWorks 教学资源）