傅里叶变换应用的例子

1. 简介

傅里叶变换就是将复杂时序函数的变为多个简单的Asin(ωx+φ)，（A为幅值，ω为频率，φ为相位）通俗易懂的话就是把复杂的化为简单。逆变换就是知道多个简单（各个频率下的相位与幅值），叠加回去就能将时序函数还原。

2. 应用

复杂的函数简称为时域函数，简单的被称为频域函数。

1. )声音的降噪：

声波是是有关于时间的时域函数，将其按频域划分，可以得到幅值与频率图。频率较低时，幅值较高的为男低音。频率中等时，幅值较高的记录了女高音的声音。当频率较高时，幅度较高处记录了噪声。此时仅需删除这些频率高的简单成分，然后通过逆变换，就能完成声音的降噪过程。

2. )图像的降噪

图像可认为是由3个通道的矩阵组成的，单通道的可以认为是黑白图。其中位置信息对应时序数据中的时间信息，图像矩阵中的每个值的大小对应时序数据中的函数值的大小。这样就能将图像信息简化为关于位置信息与该位置处的颜色深浅的值的大小的函数。当该函数采用傅里叶变换后得到频谱图是一个类似于米或十字图案。图像的轮廓部分通常是低频部分，高频往往包含了更多的细节。而图像中的噪声往往存在于高频部分，也就是说为了得到图像更清晰的轮廓，可以对低频加强，高频过滤。

3. 值得注意的是

你获得了一段时序数据，也采用fft函数获得了变换后的幅值与相位，但横坐标的频率如何获取？

这里举一个例子：（引用了Xav Zewen的文章：如何 FFT(快速傅里叶变换) 求幅度、频率（超详细 含推导过程））

你有一个最高频率f = 32kHz的模拟信号，采样频率 64kHz，对这个信号做一个16个点的FFT分析，采样点下标 n 的范围是0, 1, 2, 3, …, 15。那么64kHz的模拟频率被分成了16份，每一份是4kHz，这个4kHz被称为频率分辨率。
所以，频率图的横坐标中：
n=1 对应的f是4kHz
n=2 对应的f是8kHz
…
n=15 对应的f是60kHz
而频谱是关于n=8对称的，只需关心n = 0 ~ 7的频谱就足够了。因为，原信号的最高频率是32kHz。

频率计算公式：

f(x) = x * (Fs / n)

第0个点的频率 f(0) = 0 * (1400 / 1400) = 0
第1个点的频率 f(0) = 1 * (1400 / 1400) = 1
第2个点的频率 f(0) = 2 * (1400 / 1400) = 2
…
第200个点的频率 f(200) = 200 * (1400 / 1400) = 200
…
第1400个点的频率 f(200) = 1400 * (1400 / 1400) = 1400
（这里由于设置得很巧合，第x个点对应的频率恰好就是x）

3. Matlab代码

例子1.(完整的傅里叶变换)【周期函数】

%% 傅里叶变换
%创建一个由频率为 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) + cos(2*pi*40*t + pi/2);

% 时域信号
plot(t,x)
xlabel('时间')
ylabel('y')

% 计算信号的傅里叶变换。将变换幅值绘制为频率函数。
y = fft(x);
z = fftshift(y);
ly = length(y);
f = (-ly/2:ly/2-1)/ly*Fs;
stem(f,abs(z))
title("Double-Sided Amplitude Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("|y|")
grid

% 计算变换的相位，删除小幅值变换值。将相位绘制为频率函数。
tol = 1e-6;
z(abs(z) < tol) = 0;
theta = angle(z);
stem(f,theta/pi)
title("Phase Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("Phase/\pi")
grid

时域图

幅值频率图

相位频率图

例子2.(傅里叶变换的一半)【非周期函数】

%% FFT插值 
% 通过填充零来对信号的傅里叶变换进行插值。
% 指定信号的参数，采样频率为 80 Hz，信号持续时间为 0.8 秒。
Fs = 80;
T = 1/Fs;
L = 65;
t = (0:L-1)*T;
% 创建一个 2 Hz 正弦信号及其高次谐波的叠加。该信号包含一个 2 Hz 余弦波、一个 4 Hz 余弦波和一个 6 Hz 正弦波。
X = 3*cos(2*pi*2*t) + 2*cos(2*pi*4*t) + sin(2*pi*6*t);
% 在时域中绘制该信号。
plot(t,X)
title("Signal superposition in time domain")
xlabel("t (ms)")
ylabel("X(t)")


% 傅里叶变换
Y = fft(X);
% 计算信号的单侧幅值与频谱
f = Fs*(0:(L-1)/2)/L;
P2 = abs(Y/L);
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);

% 绘制频谱图
plot(f,P1,"-o") 
title("Single-Sided Spectrum of Original Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

时域图

幅值频率图

例子3.频率精度的提升(傅里叶变换的一半)【非周期函数】

%% 自己的函数
红外相机采样频率为30Hz,记录时长为10s，x为记录的某点的温度沿时间的函数

% 采样频率
Fs = 30;
% 采样时长10s
t = 0:1/Fs:10-1/Fs;
% 采样点的个数
L = length(t); 
% 记录时域变化
x = cos(2*pi*15*t - pi/4) + cos(2*pi*5*t + pi/2);
% 时域信号
figure
subplot(2,1,1)
plot(t,x)
xlabel('时间')
ylabel('y')

% 傅里叶变换
Y = fft(X);
% 计算信号的单侧幅值与频谱
f = Fs*(0:(L-1)/2)/L;
P2 = abs(Y/L);
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);

% 绘制频谱图
subplot(2,1,2)
plot(f,P1,"-o") 
title("Single-Sided Spectrum of Original Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

%% （增加频域的精度的操作）将尾部添0
% 为了更好地评估峰值频率，您可以通过用零填充原始信号来增加分析窗的长度。这种方法以更精确的频率分辨率
% 自动对信号的傅里叶变换进行插值。从原始信号长度确定是下一个 2 次幂的新输入长度。
% 用尾随零填充信号 X 以扩展其长度。计算填零后的信号的傅里叶变换。
n = 2^nextpow2(L);
Y = fft(x,n);
% 计算填零后的信号的单侧幅值频谱。由于信号长度 n 从 300 增加到 512，频率分辨率变为 Fs/n，即 0.5859 Hz。
f = Fs*(0:(n/2))/n;
P2 = abs(Y/L);
P1 = P2(1:n/2+1);
P1(2:end-1) = 2*P1(2:end-1);
% 绘制填零后的信号的单侧频谱。此新频谱在 0.5859 Hz 的频率分辨率内显示峰值频率。
plot(f,P1,"-o") 
title("Single-Sided Spectrum of Padded Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

图1为红外相机记录的温度与时间的数据，图2为未增加频域精度的振幅与频率图，图3为增加了频域精度的振幅与频率图。