傅里叶变换的科学应用

最新推荐文章于 2024-06-27 11:02:42 发布

大作家佚名

最新推荐文章于 2024-06-27 11:02:42 发布

阅读量1.6w

点赞数 14

分类专栏：算法/方法文章标签： fft 离散傅里叶变换功率谱功率密度滤波

本文链接：https://blog.csdn.net/wokaowokaowokao12345/article/details/69180489

版权

算法/方法专栏收录该内容

37 篇文章 6 订阅

订阅专栏

前言

最近一周一直在处理与傅里叶变换有关的事情，花了不少时间，稍微有些收获，以博文的形式记录下来希望能够自己及他人提供少许帮助。

简介

傅里叶变换可以将一个时间域信号变换到频率域。有些信号在时间域上很难看出特征，但是当转换到频率域一些特征就会变的明显。

想必当学习傅里叶变换时大家都已经知道一个例子：无限正弦波叠加就会产生一个标准的矩形波（方波或脉冲波）。如果这很难理解，就努力回忆一下大学本科高数教材中（高等数学，下册，第12章，第7节，同济版）关于傅里叶级数那一节内容。若是看到这个矩形波我们肯定想知道，组成这个矩形的那些正弦波有什么特征呢（频率、振幅、相位等）。那么傅里叶变换就可以帮我解决这个问题。

这里可能依然觉得无限正弦波叠加产生矩形波的例子有些晦涩难懂，但还有一个例子可以说是恰到好处。你一定还记得棱镜能够将白光分为红、橙、黄、绿、青、蓝、紫的例子吧。这里就可以把白光比作无限正弦波叠加后的结果，棱镜好比傅里叶变换，那么七色光就是傅里叶变换的结果。是不是经过傅里叶变换我们能知道更多的信息呢。

从书本中我们看到的最多的应该是连续函数的傅里叶变换的公式，但连续函数的傅里叶变化在计算机世界中是难以实现的。离散傅里叶变换才是我们的目标，毕竟我们都是冲着实用傅里叶变换去的，只有离散的才是眼下我们解决问题需要的。

提示：在学习过程中，有一点要注意，那就是信号的概念，切记不要狭隘的理解信号就是电磁波、声波等。一切具有波动性（非物理学中的波粒二象性）的都可以理解为信号，比如说：某条河流时间序列下的流量、某条测线一定采样间隔的高低起伏。

FFT变换计算过程

剩下时间直接进入干货阶段，本文的公式很好理解，也都是实用性公式。
http://mathworld.wolfram.com/DiscreteFourierTransform.html

离散傅里叶变换公式：

$X_k=\sum_{n=1}^{N-1}x_nW^{nk}$

傅里叶逆变换公式：

$x_k=\frac{1}{N}\sum_{n=1}^{N-1}X_nW^{-nk}$

其中， $W_n=e^{-\frac{2\pi i}{n}}$
欧拉公式：

$\begin{cases} e^{it}=cos(t)+isin(t)\\ e^{-it}=cos(t)-isin(t) \end{cases}$

实例

假设有一个序列 $x(n) = [1 2 -1 3]，n = 0, 1, 2, 3$ 。
（1）由欧拉公式， $N$ = 4 得到：

$W=e^{-2\pi/4}=cos(\frac{\pi}{2})-isin(\frac{\pi}{2})=-i$

（2）正变换结果：
$\begin{bmatrix} X(0)\\X(1)\\X(2)\\X(3) \end{bmatrix}= \begin{bmatrix} W^0&W^0&W^0&W^0\\W^0&W^1&W^2&W^3\\W^0&W^2&W^4&W^5\\W^0&W^3&W^6&W^9 \end{bmatrix} \begin{bmatrix} x(0)\\x(1)\\x(2)\\x(3) \end{bmatrix}= \begin{bmatrix} 1&1&1&1\\1&-i&-1&i\\1&-1&1&-1\\1&i&-1&-i \end{bmatrix} \begin{bmatrix} 1\\2\\-1\\3 \end{bmatrix}= \begin{bmatrix} 5\\2+i\\-5\\2-i \end{bmatrix}$

（3）逆变换结果：
$\begin{bmatrix} x(0)\\x(1)\\x(2)\\x(3) \end{bmatrix}=\frac{1}{4} \begin{bmatrix} W^0&W^0&W^0&W^0\\W^0&W^-1&W^-2&W^-3\\W^0&W^-2&W^-4&W^-5\\W^0&W^-3&W^-6&W^-9 \end{bmatrix} \begin{bmatrix} X(0)\\X(1)\\X(2)\\X(3) \end{bmatrix}=\frac{1}{4} \begin{bmatrix} 1&1&1&1\\1&i&-1&-i\\1&-1&1&-1\\1&1&-i&i \end{bmatrix} \begin{bmatrix} 5\\2+i\\-5\\2-i \end{bmatrix}= \begin{bmatrix} 1\\2\\-1\\3 \end{bmatrix}$
（4）验证，与Matlab中计算结果对比验证。
Matlab代码

clc
clear all
x = [1 2 -1 3]
y = fft(x)
x = ifft(y)

Matlab实验结果
$y = [ 5 + 0i \quad 2 + 1i \quad-5 + 0i \quad2 - 1i ]$

自己计算的结果与Matlab计算结果一致。

傅里叶变换结果意义

傅里叶变换的一个主要应用是信号分析，分析信号的基本计算包括：将双边功率谱转换为单边功率谱、调整频率精度并绘制频谱、使用 $FFT$ ，以及将功率和振幅转换为对数单位。在继续讲解前，有必要理解几个重要概念：振幅谱、功率谱、双边、单边功率谱、直流分量（DC）、对数单位。

振幅谱可由下式计算得到：

A m p l i t u d e s p e c t r u m = 2 \sqrt M a g n i t u d e [ F F T ( A ) ] N f o r i = 1 t o N 2 - 1

$Amplitude\quad spectrum= \sqrt2 \frac{Magnitude[FFT(A)]}{N} for \quad i=1\quad to\frac{N}{2}-1$

= M a g n i t u d e [ F F T ( A ) ] N f o r i = 0 (D C)

$=\frac{Magnitude[FFT(A)]}{N} for \quad i=0 \quad(DC)$
功率：可以通过平方该频率分量的幅度来获得由DFT或FFT表示的每个频率分量中的功率，既功率是振幅的平方。

功率谱（功率谱密度)：信号的功率谱密度 $powerspectraldensity(PSD)$ 是描述信号中的每单位频率对应功率的函数，功率谱密度常用单位为 $W/Hz$ 。功率谱可由下式计算：

$power spectrum \quad S_{AA}=\frac{FFT(A)FFT*(A)}{N^2}$
其中， $FFT*(A)$ 表示 $FFT(A)$ 的共轭复数。

双边、单边功率谱：傅里叶变换结果是关于直流分量对称的，既频率有正负之分，多数情况下无需再显示负极的频率信息。若使用Matlab做FFT变换，结果需要做fftshift后才是关于DC中心对称的。下文计算的功率谱是单边功率谱。

对数单位：振幅或功率谱以对数单位分贝 $(dB)$ 的形式显示。该测量单位有助于查看宽动态范围，即可在存在较大信号分量时方便地查看小信号分量。分贝是比例单位，其计算方式如下：

$dB=10\log_{10}P/P_r$
其中， $P$ 是测量功率， $P_r$ 是参考功率。使用振幅值计算分贝值：

$dB=20\log_{10}A/A_r$
其中， $A$ 是测量振幅， $A_r$ 是参考振幅。
使用振幅或功率作为同一信号的振幅平方时，结果分贝水平是完全一致的。将分贝比乘以2，等同于将比例平方。因此，无论使用振幅或功率谱，都将得到相同的分贝水平和显示。

在双边频谱中，一半的能量显示在正频率，另一半能量显示在负频率。因此，若需将双边频谱转换为单边频谱，只要舍弃数组的第二部分，并将除DC外的每个点乘以2。

直流分量：direct current component or zero-frequency component，即信号在时间域的均值。（可验证加深理解)。

原始信号有 $N$ 个采样点，经过傅里叶变换后，得到 $N$ 个点的傅里叶变换结果，即 $N$ 个复数点，每一个点对应着一个频率点。
假设采样频率为 $Fs=\frac{1}{\Delta}$ ， $\Delta$ 是采样间距。FFT之后某点 $n$ 结果用复数 $a+bi$ 表示。经过傅里叶变换，我们是要探求原始信号的一些特征（振幅、相位、功率谱等信息）。

$FFT$ 结果（除了第一个点）的振幅为其模值的 $\frac{2}{N}$ 倍。而第一个点就是直流分量，它的模值就是直流分量的 $\frac{1}{N}$ 倍。

FFT之后某点 $n$ 的模值： $A_n=\sqrt{a^2+b^2}$
频率： $F_n=(n-1)Fs/N$
振幅： $2A_n/N$ ， $n ≠1$ ，且 $n<=N/2$ ，当 $n = 1$ 时，振幅为 $A_1/N$
相位： $P_n=atan2(b,a)$ ， $atan2$ 为c/c++中求反正切的函数，其值在[-180,180]。
功率谱：对原始信号做傅里叶变换后振幅的平方。
功率谱定义：对于给定的信号，功率谱给出了落在给定频率范围内的信号功率（单位时间能量）图。

提示：功率谱，也叫做功率密度。有时候这两个概念貌似不同（不同领域），实质应该一致。出自《Numerical Recipes》数值方法（经典之作），第13.4节，第二段首句：
“The power spectrum (also called a power spectral density or PSD) .”

实例

定义信号 $S = 0.3 + 0.7*sin(2*pi*50*t) + sin(2*pi*120*t)$
从信号定义可知：直流分量为0.3（频率为0处），频率为50 $Hz$ 处振幅为0.7，
频率为120 $Hz$ 处振幅为1。离散数据特征，采样频率 $Fs$ =500，采样周期 $T=1/Fs$ ，
样本数 $N=1000$ ，时间向量： $t=(0:N-1)*T$ 。

Fs = 1000;            % 采样频率                   
T = 1/Fs;             % 采样周期       
L = 1500;             % 信号长度
t = (0:L-1)*T;        % 时间向量

S = 0.3 + 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S;
plot(t,X)
xlabel('t (seconds)')
ylabel('X(t)')
title('Signal')
figure
plot(t(1:100)*1000,X(1:100));
xlabel('t (milliseconds)')
ylabel('X(t)')

Y = fft(X);
P2 = abs(Y);
P1 = P2/(L/2);
P1(1) = P2(1)/L;
f = Fs*(0:L-1)/L;
f = f-Fs/2;
P1 = fftshift(P1);

figure
plot(f,P1) 
xlabel('f (Hz)')
ylabel('|P1(f)|')

图1为信号，图2为时间0至200 $ms$ 的信号（原信号局部图），图3为经过傅里叶变换后的频谱，从频谱中可以看到，频率为0处( $n=1$ ，数据索引)，直流分量为0.3，频率在50 $Hz$ 处( $n=101$ )振幅为0.7，频率为 $120Hz$ 处( $n=241$ )振幅为1。

频率滤波器——信号光滑处理

傅立叶滤波器是基于对信号的特定频率的一种滤波函数。它通过采用信号的傅立叶变换，然后衰减或放大特定频率，最后逆变换结果。频率滤波器的过程就是将信号从时间域经傅里叶变换转换到频率域，乘以滤波器函数，再做逆傅里叶变换转换到[空间域]。有三种基础的滤波器：低通滤波器、高通滤波器和带通滤波器。低通滤波器衰减高频，保持低频不变。空间域中的结果与平滑滤波器的结果相同。另一方面，高通滤波器在空间域中产生边缘增强或边缘检测，因为边缘包含许多高频。带通衰减极低和极高频率，保留了中频范围的频率。
$G(k,l)=F(k,l)H(k,l)$
其中 $F(k,l)$ 输入数据的傅里叶变换， $H(k,l)$ 是滤波器， $G(k,l)$ 是滤波后的数据，要获得空间域的数据还需要对其做逆傅里叶变换。
理想滤波器的反馈：

最简单的低通滤波器是理想的低通滤波器，下面滤波器看着是不是很眼熟，这就是方波的表达式。

H(k,l)=1ifk2+l2<0‾‾‾‾‾‾‾‾‾‾‾√ $H(k,l)=1\qquad if \sqrt{k^2+l^2<0}$
滤波过程可以参考下图，右下角红色部分为滤除的高拼波。傅里叶变换后信号相比原始信号变得平滑。

这里写图片描述

[低通滤波实例]

(http://stackoverflow.com/questions/28814547/matlab-low-pass-filter-using-fft)
从下图可以看到，低通滤波器在信号边缘会具有明显的失真。

clc
clear all
close all

% make our noisy function
t = linspace(1,10,1024);
x = -(t-5).^2  + 2;
y = awgn(x,0.5); 
Y = fft(y,1024);

r = 20; % range of frequencies we want to preserve

rectangle = zeros(size(Y));
rectangle(1:r+1) = 1;               % preserve low +ve frequencies
y_half = ifft(Y.*rectangle,1024);   % +ve low-pass filtered signal
rectangle(end-r+1:end) = 1;         % preserve low -ve frequencies
y_rect = ifft(Y.*rectangle,1024);   % full low-pass filtered signal

% hold on;
% plot(t,y,'g--'); 
% plot(t,x,'k','LineWidth',2); 
% plot(t,y_half,'b','LineWidth',2); 
% plot(t,y_rect,'r','LineWidth',2);
% legend('noisy signal','true signal','+ve low-pass','full low-pass','Location','southwest')

gauss = zeros(size(Y));
sigma = 8;                           % just a guess for a range of ~20
gauss(1:r+1) = exp(-(1:r+1).^ 2 / (2 * sigma ^ 2));  % +ve frequencies
gauss(end-r+1:end) = fliplr(gauss(2:r+1));           % -ve frequencies
y_gauss = ifft(Y.*gauss,1024);

hold on;
plot(t,x,'k','LineWidth',2); 
plot(t,y_rect,'r','LineWidth',2); 
plot(t,y_gauss,'c','LineWidth',2);
legend('true signal','full low-pass','gaussian','Location','southwest')

参考：
http://cn.mathworks.com/help/matlab/ref/fft.html?searchHighlight=fft&s_tid=doc_srchtitle
http://vtkvc.blog.51cto.com/1533592/314665
http://blog.csdn.net/linmingan/article/details/51077437
http://mathworld.wolfram.com/PowerSpectrum.html
https://www.wavemetrics.com/products/igorpro/dataanalysis/signalprocessing/powerspectra.htm
http://homepages.inf.ed.ac.uk/rbf/HIPR2/freqfilt.htm
http://blog.csdn.net/samkieth/article/details/49561539
http://195.134.76.37/applets/AppletFourAnal/Appl_FourAnal2.html

大作家佚名

关注

14
点赞
踩
50

收藏

觉得还不错? 一键收藏
1
评论
傅里叶变换的科学应用

前言最近一周一直在处理与傅里叶变换有关的事情，花了不少时间，稍微有些收获，以博文的形式记录下来希望能够自己及他人提供少许帮助。简介傅里叶变换可以将一个时间域信号变换到频率域。有些信号在时间域上很难看出特征，但是当转换到频率域一些特征就会变的明显。想必当学习傅里叶变换时大家都已经知道一个例子：无限正弦波叠加就会产生一个标准的矩形波（方波或脉冲波）。如果这很难理解，就努力回忆一下大学本科高数教材中（高等
复制链接

扫一扫

专栏目录