折腾来折腾去

pikipity的blog

在 Matlab 中如何做傅里叶变换

翻译自Matlab Geeks“How to Do a Fourier Transform in Matlab”

傅里叶变换是一个在很多科学和工程领域都非常有用的数学工具。傅里叶变换在信号处理、物理、通信、地质学、天文学、光学等很多领域都有应用。这个技术将一个函数或是一组数据从时域或是取样域变换到频域。这意味着,傅里叶变换可以展示一组时间序列数据的频率分量。离散傅里叶变换是将取样域的离散数据转化到频域。快速傅里叶变换是一种高效进行离散傅里叶变换的方法,并且存在很多种方法来完成快速傅里叶变换。Matlab 使用快速傅里叶变换来得到离散数据的频域分量

下面是一个在 Matlab 中如何用快速傅里叶变换来分析音频文件的例子。这个例子中的文件是记录在 note A4 (怎么翻译?不会啊)上的音叉录音。这个展示了傅里叶变换如何进行和如何在 Matlab 中使用这项技术。

%Fourier Transform of Sound File

%Load File
file = 'C:\MATLAB7\work\tuning_fork_A4';
[y,Fs,bits] = wavread(file);

Nsamps = length(y);
t = (1/Fs)*(1:Nsamps)          %Prepare time data for plot

%Do Fourier Transform
y_fft = abs(fft(y));            %Retain Magnitude
y_fft = y_fft(1:Nsamps/2);      %Discard Half of Points
f = Fs*(0:Nsamps/2-1)/Nsamps;   %Prepare freq data for plot

%Plot Sound File in Time Domain
figure
plot(t, y)
xlabel('Time (s)')
ylabel('Amplitude')
title('Tuning Fork A4 in Time Domain')

%Plot Sound File in Frequency Domain
figure
plot(f, y_fft)
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('Frequency Response of Tuning Fork A4')

音频文件 “tuning_fork_A4” 用 wavread 函数打开,这个函数返回文件的采样数据 “y"、采样频率 "Fs” 和 A/D 转换器使用的比特数 “bits"。注意,文件的扩展名 ”.wav" 不需要在使用函数的时候特别指定。就像下面将会展示的,采样频率在重建数据的时候非常重要。

快速傅里叶变换用 “fft” 函数执行。Matlab 没有 “dft” 函数因为快速傅里叶变换实际上就是计算的离散傅里叶变换。尽管快速傅里叶变换的角度在很多应用中非常中央,但是只有快速傅里叶变换的大小被保存了。“fft” 函数允许指定一个快速傅里叶变换的输出点数,但是在这个例子中,我们使用与输入点数一样数目的输出点数。在下一行,快速傅里叶变换的一半数据被舍弃了。为了这个例子的目的所以这样做了,但是在很多应用中,整个波谱都会用到(译者注:我认为这里舍弃一半的点是因为FFT是关于采样频率的一半对称的,所以只要看一半就可以了)。接下来的一行,将会用于横坐标的数据通过使用采样频率和时遇采样数量准备好了。这一步对于确定包含在音频文件的实际频率是很重要的。

接下来,原始数据在时域上被画了下来,快速傅里叶变换的数据也被画了出来。为了展示在峰值的频率上的更多详情,在这个画中,x轴被限定在了 [0,1000] 的范围中。注意,在大约440Hz处,频率响应有一个峰值,这个就是 “note A4” 的频率。在其他频率也有一些很小东西,这些估计是音叉了。对于其他乐器比如吉他,在频率响应的峰值的整数倍上都可以看见谐波。

傅里叶变换在很多不同的领域都是很有用的工具。二维的傅里叶变换也常常用在图像上。尝试一下上面的代码,看看你能不能得到一样的结果。

PS:有人问 FFT 结果的幅值问题。本质来说,FFT结果的幅值单位是什么并不重要,只要你在分析过程中,需要分析的两个幅值的单位是统一的就可以了。Matlab 的 FFT 最终结果的绝对值的确看上去并不好看(太大了),根据 Matlab 帮助文件,FFT 的最终结果还需要进行一下小调整再来使用比较好(上面的图未作调整),这里有我根据 Matlab 帮助文档编写的一个计算 FFT 并绘制 FFT 结果幅值图的函数,具体如何使用,请看函数内的注释说明。



Comments