前言
最近一周一直在处理与傅里叶变换有关的事情,花了不少时间,稍微有些收获,以博文的形式记录下来希望能够自己及他人提供少许帮助。
简介
傅里叶变换可以将一个时间域信号变换到频率域。有些信号在时间域上很难看出特征,但是当转换到频率域一些特征就会变的明显。
想必当学习傅里叶变换时大家都已经知道一个例子:无限正弦波叠加就会产生一个标准的矩形波(方波或脉冲波)。如果这很难理解,就努力回忆一下大学本科高数教材中(高等数学,下册,第12章,第7节,同济版)关于傅里叶级数那一节内容。若是看到这个矩形波我们肯定想知道,组成这个矩形的那些正弦波有什么特征呢(频率、振幅、相位等)。那么傅里叶变换就可以帮我解决这个问题。
这里可能依然觉得无限正弦波叠加产生矩形波的例子有些晦涩难懂,但还有一个例子可以说是恰到好处。你一定还记得棱镜能够将白光分为红、橙、黄、绿、青、蓝、紫的例子吧。这里就可以把白光比作无限正弦波叠加后的结果,棱镜好比傅里叶变换,那么七色光就是傅里叶变换的结果。是不是经过傅里叶变换我们能知道更多的信息呢。
从书本中我们看到的最多的应该是连续函数的傅里叶变换的公式,但连续函数的傅里叶变化在计算机世界中是难以实现的。离散傅里叶变换才是我们的目标,毕竟我们都是冲着实用傅里叶变换去的,只有离散的才是眼下我们解决问题需要的。
提示:在学习过程中,有一点要注意,那就是信号的概念,切记不要狭隘的理解信号就是电磁波、声波等。一切具有波动性(非物理学中的波粒二象性)的都可以理解为信号,比如说:某条河流时间序列下的流量、某条测线一定采样间隔的高低起伏。
FFT变换计算过程
剩下时间直接进入干货阶段,本文的公式很好理解,也都是实用性公式。
http://mathworld.wolfram.com/DiscreteFourierTransform.html
离散傅里叶变换公式:
Xk=∑N−1n=1xnWnk
傅里叶逆变换公式:
xk=1N∑N−1n=1XnW−nk
其中,
Wn=e−2πin
欧拉公式:
{eit=cos(t)+isin(t)e−it=cos(t)−isin(t)
实例
假设有一个序列
x(n)=[12−13],n=0,1,2,3
。
(1)由欧拉公式,
N
= 4 得到:
(2) 正变换结果:
⎡⎣⎢⎢⎢⎢X(0)X(1)X(2)X(3)⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢W0W0W0W0W0W1W2W3W0W2W4W6W0W3W5W9⎤⎦⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢x(0)x(1)x(2)x(3)⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢11111−i−1i1−11−11i−1−i⎤⎦⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢12−13⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢52+i−52−i⎤⎦⎥⎥⎥⎥
(3)逆变换结果:
⎡⎣⎢⎢⎢⎢x(0)x(1)x(2)x(3)⎤⎦⎥⎥⎥⎥=14⎡⎣⎢⎢⎢⎢W0W0W0W0W0W−1W−2W−3W0W−2W−4W−6W0W−3W−5W−9⎤⎦⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢X(0)X(1)X(2)X(3)⎤⎦⎥⎥⎥⎥=14⎡⎣⎢⎢⎢⎢11111i−111−11−i1−i−1i⎤⎦⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢52+i−52−i⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢12−13⎤⎦⎥⎥⎥⎥
(4)验证,与Matlab中计算结果对比验证。
Matlab代码
clc
clear all
x = [1 2 -1 3]
y = fft(x)
x = ifft(y)
Matlab实验结果
y=[5+0i2+1i−5+0i2−1i]
自己计算的结果与Matlab计算结果一致。
傅里叶变换结果意义
傅里叶变换的一个主要应用是信号分析,分析信号的基本计算包括:将双边功率谱转换为单边功率谱、调整频率精度并绘制频谱、使用 FFT ,以及将功率和振幅转换为对数单位。在继续讲解前,有必要理解几个重要概念:振幅谱、功率谱、双边、单边功率谱、直流分量(DC)、对数单位。
振幅谱可由下式计算得到:
功率:可以通过平方该频率分量的幅度来获得由DFT或FFT表示的每个频率分量中的功率,既功率是振幅的平方。
功率谱(功率谱密度):信号的功率谱密度 powerspectraldensity(PSD) 是描述信号中的每单位频率对应功率的函数,功率谱密度常用单位为 W/Hz 。功率谱可由下式计算:
powerspectrumSAA=FFT(A)FFT∗(A)N2
其中,
FFT∗(A)
表示
FFT(A)
的共轭复数。
双边、单边功率谱:傅里叶变换结果是关于直流分量对称的,既频率有正负之分,多数情况下无需再显示负极的频率信息。若使用Matlab做FFT变换,结果需要做fftshift后才是关于DC中心对称的。下文计算的功率谱是单边功率谱。
对数单位:振幅或功率谱以对数单位分贝 (dB) 的形式显示。该测量单位有助于查看宽动态范围,即可在存在较大信号分量时方便地查看小信号分量。分贝是比例单位,其计算方式如下:
dB=10log10P/Pr
其中,
P
是测量功率,
dB=20log10A/Ar
其中,
A
是测量振幅,
使用振幅或功率作为同一信号的振幅平方时,结果分贝水平是完全一致的。将分贝比乘以2,等同于将比例平方。因此,无论使用振幅或功率谱,都将得到相同的分贝水平和显示。
在双边频谱中,一半的能量显示在正频率,另一半能量显示在负频率。因此,若需将双边频谱转换为单边频谱,只要舍弃数组的第二部分,并将除DC外的每个点乘以2。
直流分量:direct current component or zero-frequency component,即信号在时间域的均值。(可验证加深理解)。
原始信号有
N
个采样点,经过傅里叶变换后,得到
假设采样频率为
FFT 结果(除了第一个点)的振幅为其模值的 2N 倍。而第一个点就是直流分量,它的模值就是直流分量的 1N 倍。
FFT之后某点
n
的模值:
频率:
Fn=(n−1)Fs/N
振幅:
2An/N
,
n≠1
,且
n<=N/2
,当
n=1
时, 振幅为
A1/N
相位:
Pn=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‾‾‾‾‾‾‾‾‾‾‾√
滤波过程可以参考下图,右下角红色部分为滤除的高拼波。傅里叶变换后信号相比原始信号变得平滑。
[低通滤波实例]
(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