网站开发人员需求,平面设计赚钱平台,乐清网站改版公司,wordpress设主题前言
最近复现音乐驱动舞蹈的文章《Dancing-to-Music Character Animation》#xff0c;用到了与傅里叶变换很相似的称为常Q变换的方法去分割音乐#xff0c;所以对傅里叶变换做了一个小了解#xff0c;本文不深入各种乱糟糟的理论#xff0c;比如什么蝶形算法啥的#x…前言
最近复现音乐驱动舞蹈的文章《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
#简介
傅里叶变换意思就是任何一个输入信号都可以使用多个余弦波叠加而成简单的说就是把时序信号转换成频域信息。我们最终需要的就是找到这些余弦波的相关参数幅值相位。
复习一下三角函数的标准式 yAcos(wxb)kyAcos(wxb)k yAcos(wxb)k AAA代表振幅函数周期是2πw\frac{2\pi}{w}w2π频率是周期的倒数w2π\frac{w}{2\pi}2πwbbb是函数初相位kkk在信号处理中称为直流分量。
在很多工具的实现中余弦波的个数就是信号长度但是在理论公式中会出现一个参数N是采样点通常称为N点FFT。 X(k)∑n1Nx(n)∗exp(−−1∗2π∗(k−1)∗(n−1)/N),1kN.X(k) \sum _{n1}^N x(n)*\exp(-\sqrt{-1}*2\pi*(k-1)*(n-1)/N), 1 k N. X(k)n1∑Nx(n)∗exp(−−1∗2π∗(k−1)∗(n−1)/N),1kN. 上述公式就是DFT的求解方法了不考虑它的优化方法FFT我们直接在matlab中码公式最后与FFT的结果做对比即可验证写的算法对不对。
#实例
假设我们的输入信号是函数是 S0.20.7∗cos(2π∗50t20180π)0.2∗cos(2π∗100t70180π)S0.20.7*\cos(2\pi *50t\frac{20}{180}\pi) 0.2*\cos(2\pi*100t\frac{70}{180}\pi) S0.20.7∗cos(2π∗50t18020π)0.2∗cos(2π∗100t18070π) 可以发现直流分量是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.20.7*cos(2*pi*50*t20/180*pi) 0.2*cos(2*pi*100*t70/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/21);
P1(2:end-1) 2*P1(2:end-1);首先直接对原始信号进行傅里叶变换得到YYY它包含实部与虚部然后求解归一化后YYY各项的模得到P2P2P2由于matlab求解的结果包含对称的两个频谱称为双侧频谱我们只需要取一半得到P1P1P1此时只需要将除第一个元素外的元素乘以2即可得到幅值信息其实求解的根本就是来源于YYYYYY有多少项就说明求解了多少个叠加的余弦波接下来解释如何求解各参数
直流分量就是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的位置分别是050100附近其实我们去看具体的位置发现是1,51,101此三个位置的YYY值分别为 Y(1)ans 200.0000 Y(51)ans 3.2889e02 1.1971e02i Y(101)ans 34.2020 93.9693i那么按照描述我们得到 直流分量Y(1)Fs0.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−π)也就是说如果幅值是负数那么只需要在正数幅值的基础上变化一初相位。比如把例子的函数换成 S0.2−0.7∗cos(2π∗50t20180π)0.2∗cos(2π∗100t70180π)S0.2-0.7*\cos(2\pi *50t\frac{20}{180}\pi) 0.2*\cos(2\pi*100t\frac{70}{180}\pi) S0.2−0.7∗cos(2π∗50t18020π)0.2∗cos(2π∗100t18070π) 那么求解频率50对应余弦波的赋值为0.7初相位为-160
IFFT变换
顾名思义IFFT就是逆傅里叶变换用Y重构信号其实我们通过Y都已经分析出了原始信号所需的余弦波的各参数肯定能重构原始数据这里还是做个实验吧用IFFT函数
figure
pred_Xifft(Y);
plot(pred_X,r-)
hold on
plot(S,b*)DFT变换
按照公式手撸DFT看看计算得到YYY与matlab自带的FFT结果是否一致 X(k)∑n1Nx(n)∗exp(−−1∗2π∗(k−1)∗(n−1)/N),1kN.X(k) \sum _{n1}^N x(n)*\exp(-\sqrt{-1}*2\pi*(k-1)*(n-1)/N), 1 k N. X(k)n1∑Nx(n)∗exp(−−1∗2π∗(k−1)∗(n−1)/N),1kN.
%% DFT变换
% N
% X(k) sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 k N.
% n1
DFT_Xzeros(1,L);
NL;
for k1:Lfor n1:NDFT_X(k)DFT_X(k)S(n)*exp(-1j*2*pi*(k-1)*(n-1)/N);end
end
P2abs(DFT_X/L);
P1 2*P2(1:L/21);
f Fs*(0:(L/2))/L;
figure
plot(f,P1)
xlabel(频率)
ylabel(振幅)
title(DFT变换)再计算与FFT求解的结果的误差
%% FFT变换
FFT_Xfft(S);
figure
plot(abs(FFT_X-DFT_X))
title(DFT和FFT误差)IDFT
同样按照公式手撸逆离散傅里叶变换 x(n)1N∑k1NX(k)∗exp(−1∗2π∗(k−1)∗(n−1)/N),1nN.x(n) \frac{1}{N}\sum _{k1}^N X(k)*\exp(\sqrt{-1}*2\pi*(k-1)*(n-1)/N), 1 n N. x(n)N1k1∑NX(k)∗exp(−1∗2π∗(k−1)∗(n−1)/N),1nN.
%% IDFT变换
% N
% x(n) (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 n N.
% k1
DFT_rec_xzeros(1,k);
for n1:Lfor k1:LDFT_rec_x(n)DFT_rec_x(n)DFT_X(k)*exp( 1j*2*pi*(k-1)*(n-1)/N);end
end
DFT_rec_xDFT_rec_x./N;
rec_errreal(DFT_rec_x)-S;
figure
plot(rec_err)
title(IFFT数据重构误差)与IFFT的结果对比一下
%% IFFT变换
FFT_rec_xifft(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的蝶形算法_
博客已同步更新到微信公众号中有兴趣可关注一波微信公众号与博客同步更新个人的学习笔记
