前言

最近复现音乐驱动舞蹈的文章《Dancing-to-Music Character Animation》，用到了与傅里叶变换很相似的称为常Q变换的方法去分割音乐，所以对傅里叶变换做了一个小了解，本文不深入各种乱糟糟的理论，比如什么蝶形算法啥的，单纯讨论离散傅里叶变换(DFT)，我们常见的是快速傅里叶变换(FFT)，其实FFT是对DFT的一个计算优化，主要是剔除DFT里面有周期性之类的被冗余计算，但是FFT的算法有点小复杂，有兴趣深入理论的请移步如下几篇博文：

如何理解和掌握快速傅里叶变换的计算和概念？

如何通俗地解释什么是离散傅里叶变换？

傅里叶分析之掐死教程（完整版）更新于2014.06.06

傅里叶变换

快速傅里叶变换

第三章离散傅里叶变换(DFT) 及其快速算法(FFT)

傅里叶级数和傅里叶变换是什么关系？

如何理解傅里叶变换公式？

An Introduction to the Fast Fourier Transform

FFT的详细解释，相信你看了就明白了

如果想仔细学习FFT，最好是看完上述推荐的博文并额外找资料，我是不想看了，因为我看着看着发现自己掉头发了o(╯□╰)o

#简介

傅里叶变换意思就是任何一个输入信号都可以使用多个余弦波叠加而成，简单的说就是把时序信号转换成频域信息。我们最终需要的就是找到这些余弦波的相关参数：幅值，相位。

复习一下三角函数的标准式：
y=Acos(wx+b)+ky=Acos(wx+b)+k y=Acos(wx+b)+k
AAA代表振幅，函数周期是2πw\frac{2\pi}{w}w2π，频率是周期的倒数w2π\frac{w}{2\pi}2πw，bbb是函数初相位，kkk在信号处理中称为直流分量。

在很多工具的实现中，余弦波的个数就是信号长度，但是在理论公式中会出现一个参数N，是采样点，通常称为N点FFT。
X(k)=∑n=1Nx(n)∗exp⁡(−−1∗2π∗(k−1)∗(n−1)/N),1<=k<=N.X(k) = \sum _{n=1}^N x(n)*\exp(-\sqrt{-1}*2\pi*(k-1)*(n-1)/N), 1 <= k <= N. X(k)=n=1∑Nx(n)∗exp(−−1∗2π∗(k−1)∗(n−1)/N),1<=k<=N.
上述公式就是DFT的求解方法了，不考虑它的优化方法FFT，我们直接在matlab中码公式，最后与FFT的结果做对比即可验证写的算法对不对。

#实例

假设我们的输入信号是函数是
S=0.2+0.7∗cos⁡(2π∗50t+20180π)+0.2∗cos⁡(2π∗100t+70180π)S=0.2+0.7*\cos(2\pi *50t+\frac{20}{180}\pi) + 0.2*\cos(2\pi*100t+\frac{70}{180}\pi) S=0.2+0.7∗cos(2π∗50t+18020π)+0.2∗cos(2π∗100t+18070π)
可以发现直流分量是0.2，以及两个余弦函数的叠加，余弦函数的幅值分别为0.7和0.2，频率分别为50和100，初相位分别为20度和70度。

画出信号图：

Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 1000;             % Length of signal
t = (0:L-1)*T;        % Time vector
S = 0.2+0.7*cos(2*pi*50*t+20/180*pi) + 0.2*cos(2*pi*100*t+70/180*pi) ;
plot(1000*t(1:50),S(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')

FFT变换

先使用matlab默认的快速傅里叶变换函数FFT求解叠加余弦的各参数

Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

首先直接对原始信号进行傅里叶变换得到YYY，它包含实部与虚部，然后求解归一化后YYY各项的模得到P2P2P2，由于matlab求解的结果包含对称的两个频谱，称为双侧频谱，我们只需要取一半得到P1P1P1，此时只需要将除第一个元素外的元素乘以2即可得到幅值信息，其实求解的根本就是来源于YYY，YYY有多少项，就说明求解了多少个叠加的余弦波，接下来解释如何求解各参数：

直流分量：就是Y的第一个值除以采样频率，一般来说是非复数
频率：采样频率采样长度∗(第几项−1)\frac{采样频率}{采样长度}*(第几项-1)采样长度采样频率∗(第几项−1)，本例中采样频率是1000，长度也是1000，那么YYY的第二项频率就是1，第三项频率是2，其实最终情况下，我们选取频率不接近0的值。
幅值：2∗abs(Y各项采样长度)2*abs(\frac{Y各项}{采样长度})2∗abs(采样长度Y各项)
初相位：atan2(Y的虚部Y的实部)atan2(\frac{Y的虚部}{Y的实部})atan2(Y的实部Y的虚部)转角度制表示

从P1P1P1的图中，我们很容易看出来幅值不接近0的位置分别是0，50，100附近，其实我们去看具体的位置发现是1,51,101，此三个位置的YYY值分别为：

>> Y(1)ans =200.0000>> Y(51)ans =3.2889e+02 + 1.1971e+02i>> Y(101)ans =34.2020 +93.9693i

那么按照描述，我们得到：

直流分量：Y(1)Fs=0.2\frac{Y(1)}{Fs}=0.2FsY(1)=0.2
频率：第51项和101项的频率分别为50和100
幅值：2∗abs(Y(51)L)=0.72*abs\left(\frac{Y(51)}{L}\right)=0.72∗abs(LY(51))=0.7 同理2∗abs(Y(101)L)=0.22*abs\left(\frac{Y(101)}{L}\right)=0.22∗abs(LY(101))=0.2

初相位：

>> rad2deg(atan2(imag(Y(51)),real(Y(51))))ans =20.0000>> rad2deg(atan2(imag(Y(101)),real(Y(101))))ans =70.0000

【注1】: 在实际应用中，一般让余弦波的的数量与信号长度相同，如果不相同，那就是理论中常说的N点FFT

【注2】：幅值是通过模求解的，但是模一般都是正数，如果叠加信号的幅值是复数呢？不要忘记了−cos(x)=cos(x−π)-cos(x)=cos(x-\pi)−cos(x)=cos(x−π)，也就是说如果幅值是负数，那么只需要在正数幅值的基础上变化一初相位。比如把例子的函数换成：
S=0.2−0.7∗cos⁡(2π∗50t+20180π)+0.2∗cos⁡(2π∗100t+70180π)S=0.2-0.7*\cos(2\pi *50t+\frac{20}{180}\pi) + 0.2*\cos(2\pi*100t+\frac{70}{180}\pi) S=0.2−0.7∗cos(2π∗50t+18020π)+0.2∗cos(2π∗100t+18070π)
那么求解频率50对应余弦波的赋值为+0.7，初相位为-160

IFFT变换

顾名思义，IFFT就是逆傅里叶变换，用Y重构信号，其实我们通过Y都已经分析出了原始信号所需的余弦波的各参数，肯定能重构原始数据，这里还是做个实验吧，用IFFT函数：

figure
pred_X=ifft(Y);
plot(pred_X,'r-')
hold on
plot(S,'b*')

DFT变换

按照公式手撸DFT，看看计算得到YYY与matlab自带的FFT结果是否一致
X(k)=∑n=1Nx(n)∗exp⁡(−−1∗2π∗(k−1)∗(n−1)/N),1<=k<=N.X(k) = \sum _{n=1}^N x(n)*\exp(-\sqrt{-1}*2\pi*(k-1)*(n-1)/N), 1 <= k <= N. X(k)=n=1∑Nx(n)∗exp(−−1∗2π∗(k−1)∗(n−1)/N),1<=k<=N.

%% DFT变换
%                    N
%      X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                   n=1
DFT_X=zeros(1,L);
N=L;
for k=1:Lfor n=1:NDFT_X(k)=DFT_X(k)+S(n)*exp(-1j*2*pi*(k-1)*(n-1)/N);end
end
P2=abs(DFT_X/L);
P1 = 2*P2(1:L/2+1);
f = Fs*(0:(L/2))/L;
figure
plot(f,P1)
xlabel('频率')
ylabel('振幅')
title('DFT变换')

再计算与FFT求解的结果的误差

%% FFT变换
FFT_X=fft(S);
figure
plot(abs(FFT_X-DFT_X))
title('DFT和FFT误差')

IDFT

同样按照公式手撸逆离散傅里叶变换
x(n)=1N∑k=1NX(k)∗exp⁡(−1∗2π∗(k−1)∗(n−1)/N),1<=n<=N.x(n) = \frac{1}{N}\sum _{k=1}^N X(k)*\exp(\sqrt{-1}*2\pi*(k-1)*(n-1)/N), 1 <= n <= N. x(n)=N1k=1∑NX(k)∗exp(−1∗2π∗(k−1)∗(n−1)/N),1<=n<=N.

%% IDFT变换
%                    N
%      x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
%                   k=1
DFT_rec_x=zeros(1,k);
for n=1:Lfor k=1:LDFT_rec_x(n)=DFT_rec_x(n)+DFT_X(k)*exp( 1j*2*pi*(k-1)*(n-1)/N);end
end
DFT_rec_x=DFT_rec_x./N;
rec_err=real(DFT_rec_x)-S;
figure
plot(rec_err)
title('IFFT数据重构误差')

与IFFT的结果对比一下：

%% IFFT变换
FFT_rec_x=ifft(FFT_X);
figure
plot(abs(DFT_rec_x-FFT_rec_x))
title('IDFT和IFFT误差')

后记

折腾了这么多，其实就是为了一个字：懒。为了避免复杂的FFT的理论理解，我直接按照DFT的公式码代码，取得了与FFT一样的结果，对于不喜欢复杂理论的同志，可以在代码实现FFT的时候直接采用DFT的代码作为替代品，虽然时间复杂度增大很多，但是理论理解起来简单很多倍不是。

贴一下文章代码，此外还有一个FFT的手动实现：链接：https://pan.baidu.com/s/1mUdclEQ3tgQUvZQ4XffN3g 密码：08tg

等我不掉头发了，再去纠结FFT的蝶形算法^_

博客已同步更新到微信公众号中，有兴趣可关注一波，微信公众号与博客同步更新个人的学习笔记

【音频处理】离散傅里叶变换相关推荐

C++实现二维离散傅里叶变换
在上一篇文章<C++实现一维离散傅里叶变换>中,我们介绍了一维信号傅立叶变换的公式和C++实现,并阐述了频域幅值的意义. 一维傅立叶变换只适用于一维信号,例如音频数据.心脑电图等. 在图像 ...
（数字图像处理MATLAB+Python）第四章图像正交变换-第一节：离散傅里叶变换
文章目录一:一维离散傅里叶变换 (1)定义 (2)实例二:一维快速傅里叶变换 (1)定义 (2)实例三:二维离散傅里叶变换 (1)定义 (2)程序四:二维离散傅里叶变换的性质 (1)可分性 ( ...
《OpenCV3编程入门》学习笔记5 Core组件进阶（五）离散傅里叶变换（DFT）
第5章 Core组件进阶 5.5 离散傅里叶变换(Discrete Fourier Transform,DFT) 5.5.1 离散傅里叶变换原理 1.对一张图像使用傅里叶变换就是把它分解成正弦和余弦, ...
独家｜OpenCV 1.7 离散傅里叶变换
翻译:陈之炎校对:李海明本文约2400字,建议阅读5分钟本文为大家介绍了OpenCV离散傅里叶变换. 目标本小节将寻求以下问题的答案: 什么是傅立叶变换,为什么要使用傅立叶变换? 如何在OpenC ...
离散傅里叶变换（DFT）（为了使用而学习的DFT）
1. 离散周期信号的傅里叶级数及其系数(DFS) 1)针对的对象:周期离散序列,设周期为N: 2)像连续周期信号那样用傅里叶级数表示信号,也即周期序列x[n]的傅里叶级数(DFS)表示: 其中: 从上 ...
用离散傅里叶变换来实现OFDM
这就是使用离散傅里叶变换来实现OFDM的原理图: 那分帧分组,编码映射是什么情况呢? 用途是将输入的比特流先分帧,然后在帧中分组,然后再串并转换: 然后再将这些分组编码,并使编码与QAM或QPSK的星 ...
ubuntu 使用FFTW快速计算离散傅里叶变换
FFTW ( the Faster Fourier Transform in the West) 是一个快速计算离散傅里叶变换的标准C语言程序集,其由MIT的M.Frigo 和S. Johnson 开 ...
Opencv 实现图像的离散傅里叶变换（DFT）、卷积运算（相关滤波）
原文:http://blog.csdn.net/ikerpeng/article/details/41845545?utm_source=tuicool&utm_medium=referral ...
基于OpenCV完成离散傅里叶变换
基于OpenCV完成离散傅里叶变换目标学会使用函数: cv::copyMakeBorder() , cv::merge() , cv::dft() , cv::getOptimalDFTSize( ...

【音频处理】离散傅里叶变换

前言

FFT变换

IFFT变换

DFT变换

IDFT

后记

【音频处理】离散傅里叶变换相关推荐

最新文章

热门文章