2020-12-13

最新推荐文章于 2022-08-17 16:21:57 发布

Darren_R

最新推荐文章于 2022-08-17 16:21:57 发布

阅读量389

点赞数 1

文章标签：傅立叶分析

本文链接：https://blog.csdn.net/alinasina/article/details/111088431

版权

一、Introduction

为什么傅里叶感觉很复杂？

一是名字就很搞人

什么连续傅里叶、离散傅里叶、短时傅里叶、快速傅里叶等等

最重要的一点，变量的不统一，比如常见的（冲激，冲击、脉冲）属于翻译搞人；再比如冲激串函数（也叫dirac comb function)，这个函数下标有用T(代表周期）、P（代表脉冲）、S(代表采样）‘

信号采样一般我们认为是得到数字，这个也叫序列，但其实得到的还有可能是函数，叫做抽样后的函数，多亏了冲激函数串的引入可以

预备概念：

n次单位根：

傅里叶各种变换首先是来源于傅里叶先生提出的那句惊世骇俗的话：任意信号均可以由不同频率的正弦波构成，当不同频率是指基波的正整数倍。那么为什么是正整数呢，正符合现实，整数符合三角波的正交性。

二、代码验证

CTFS(Continous Time Fourier Series 连续时间傅里叶级数)

单边频谱计算系数ak,bk和a0的代码(k代表第k次谐波，也可以理解为“对频率采样”，使用k是为了与n区分，n一般被理解为"对时间采样”）

下面是一个矩形信号的CTFS,注意：后面的画图我们均采用stem函数，相比于plot直观多了，常用于离散信号的画图；但要注意这个stem画出来的可不是冲激串函数，只是长得一样。

syms t k;
T=10;tao=5;E=1;
k=1:1:5;%从1次谐波求到5次谐波
y=E*(heaviside(t+tao/2)-heaviside(t-tao/2));%y表示一个周期内的信号函数，并非整个周期的信号函数

a0=1/T*int(y,t,-T/2,T/2); % a0
ak=2/T*int(y*cos(2*pi*k*t/T),t,-T/2,T/2); % ak
bk=2/T*int(y*sin(2*pi*k*t/T),t,-T/2,T/2); % bk

stem(k,ak);
hold on;
stem(0,a0);
grid on;
title('周期矩形波信号的CTFS-单边频谱');
%result:a0=1/2;ak[5]=[2/pi,0,-2/(3*pi),0,2/(5*pi)],bk[k]=[0,0,0,0,0];

但是上面这种计算CTFS-单边频谱的方法不是我们常用的，我们常计算双边频谱；而且直接对符号变量积分太硬核，不如换成数值积分，轻松省事~

一种简单的计算双边频谱的方法是利用上述代码再结合ck=1/2*(ak-jbk)来计算

另一种直接套定义去算ck

syms t k;
T=10;tao=5;E=1;
k=-5:1:5;%从-5次谐波求到5次谐波
y=E*(heaviside(t+tao/2)-heaviside(t-tao/2));%y表示一个周期内的信号函数，并非整个周期的信号函数

ck=1/T*int(y*exp(-1j*2*pi*k*t/T),t,-T/2,T/2); % ck

stem(k,ck);
grid on;
title('周期矩形波信号的CTFS-双边频谱');
%result:ck[11]=[1/(5*pi),0,-1/(3*pi),0,1/pi,1/2,1/pi,0,-1/(3*pi),0,1/(5*pi)];
t=-T/2:0.01:T/2;

其实吧ck相对于a0,ak,bk的表现形式好处很多呀，一个是直流量a0不用单独计算了，其次用复指数替代三角函数（欧拉公式）的话整个公式简洁多了’；

ck的数域范围相比于a0,ak,bk数域范围从实数拓宽到复数域，复数嘛，自然就有了幅值和相角；

由于负频率的出现（人为引入），ck多了另一半谱，右半平面谱相对于左半平面谱来说，幅值谱关于y轴对称，相位谱关于原点对称，所以吧，我们一般就看右半平面谱就好了。所以我们常说的傅里叶频谱就是一个有半平面的“复数”谱。

上面的一个周期内的矩形信号由两个阶跃函数构成，我们称这种方法是解析的方法；现在用数值的方法构建这个矩形信号

CTFS和CTFT由于时间的连续，因此工程上选取了离散时间的傅里叶变换

根据对称性，就有DTFS和DTFT，DTFS还好啦，但是DTFT问题最严重。我们知道离散信号的周期是N，代表我采样了N个点。它的谐波最高也是N。DTFT和CTFT一样，都是周期趋于无穷得到的。但是DTFT根本没法用的原因在于它要在时域采样无穷个点，而且这要做出的结果频域还仍然是连续的，果然不行。DFT干脆说，我不采无穷多个点了，就采N个，然后做DTFT，做出来的效果和N个点做周期延拓再进行DTFS效果一样。

所以就有了这两句话，DFT是DTFT在频域采样后的结果；DFT也是时域采样进行周期延拓后进行DTFS的结果。所以我们就只关注DFT或者DTFS的实现了。

但是不管怎么做，我们希望做到对任意一个信号的频谱都能从工程上实际还原回去，验证的标准可以是理论计算和实际的对比。

DTFS（DIscrete Time Fourier Series,离散时间傅里叶级数）

代码实现：

function Xk = dtfs(xn) 
%自定义函数
%脚本文件:dtfs.m
%dtfs函数：discrete time fourier series 离散时间傅里叶级数
%定义变量：
%xn---向量，表征采样序列
%Xk---向量
%版本记录：2020-12-20 Darren
N=length(xn);% 序列xn的长度
n=[0:1:N-1];% n代表N个样本中的第n个采样
k=n;% k代表第k次谐波，谐波次数最多为N-1次
WN = exp(-1i*2*pi/N);% WN是单位根的一种表示方法
nk=n'*k;% 这里代表dtfs公式中的指数
WNnk=WN.^nk;% 代表dtfs公式中第nk次的单位根
Xk =xn*WNnk;% 代表第k次谐波的大小，它是是一个(复指）数，有模和相角
end

Ex1:正弦信号取样后进行dtfs

正弦信号的生成

%正弦信号生成与测试代码，需要调用dtfs/或者fft
clear;
ts=1/32;%取样时间
fs=1/ts;%周期
N=35;%总取样次数
n=0:N-1;
xn=cos(2*pi*n*ts);%取离散信号数据
subplot(2,1,1);
stem(n,xn);%N=32,时域取样图
title(['N=',num2str(N),'时域取样图']);
xlabel('n(第n个采样)');
x1=dtfs(xn);
subplot(2,1,2);
stem(n,abs(x1));%绘制时域取样图
title('频域取样图');
xlabel('k(第k次谐波)');

按照公式来看，cos2*pi*t做CTFS应该出现的幅值为pi，正负频率为-2*pi和2*pi的基波频率，但是我们做dtfs(或dft，fft)的时候发现在非基波频率上也会出现值，而不是0，这是为什么呢？一开始我以为是频率混叠，后面推导后发现真正的原因在于采样点数N和采样时间Ts的问题。以下图为例

可以很明显的看到造成的N=32是我们理想的波形，但是当N=25或N=35时在其余次谐波上均有值。当然从时域可以看出N=25和N=35的波形采样不是一个正弦波，而且如果做周期延拓那就更不像了。所以自然频谱就是正弦波的那两条线了。【补充：后面看到一些资料提到，上面这个现象叫做“频率泄露”，就是说本来只应该只有一个频率的地方却在其附近出现了频率

细心的你一定会发现频谱为什么是不是关于k=0对称呢，反而像是关于k=某个值对称呢；关键的原因在于我们的公式，我们k的取值范围是[0，N-1]，而不是[-(N-1)/2,（N-1)/2-1]。所以我们心目中双边频谱左半平面的频谱就“右移拼接到了最右边了

我们常说的fft是dft的改良，计算乘法的次数少了。我们增加了一组用matlab自带fft算法计算的结果

这样再次证明自己编写的dft/dtfs代码和matlab给出结果是一致的。

但是我们这里还有一个问题，那就是上面每一张图的横轴要么是n，要么是k，由于dft算法本身不需要知道时间信息，所以要想知道时间就要进行一些时间上的变化

谈到FFT，我们很多人并不知道要选多个点做FFT合适，如何采样？

首先为了方便N取2的整数次方

为了得到fft后的真实幅值，需要做一个乘以2除以N

为什么fft需要加窗，频率泄露是什么鬼【26】

自己还有个小问题待解决？拉氏变换后s=jw得到的傅里叶变换到频谱（bode图）其幅值如何变换为真实的幅值

测频谱为什么要去掉直流，直流对低频的频谱有影响？

频率细分法和补0又是什么鬼？

Ex2: 方波信号的采样

方法一是利用square函数或者阶跃函数取样获得，方法2是使用ones和zeros的矩阵直接生成（但是这种方法生成方波的一个缺陷是点数这样的波形一定是一个梯形波，在间断点处有斜坡，哎，这也是现实和理想的差距）

短时傅里叶变换

附录代码：

常用的画图函数

subplot(2,2,1),plot(t1,f1),hold on;plot(t1,y,'r:');title('周期矩形波的形成—基波'),axis([-2.5,2.5,-0.5,1.1])
title('周期方波信号','Fontsize',8);
xlabel('t   (单位:s)', 'Fontsize',8);

离散数学里常用的matlab函数

F1=int(f,x,a,b) — 对函数f(x)在（a,b)区间求定积分,定积分大小是F1
F2=subs(s,OLD,NEW)－用新变量NEW代替S中的指定变量OLD。
F3=vpa(x,n) ： 显示可变精度计算；x为符号变量，n表示要精确计算的位数。

常用信号生成代码

%单位阶跃信号
t=-10:0.001:10;
y=(t>=0);
y=heaviside(t);
plot(t,y,'r');

%单位冲激信号
x=-100:0.1:100;
y=dirac(x); %狄拉克函数
y=5*sign(y); %尺度变换，否则显示不出 infinity
plot(x,y,'r')

%指数信号
A=1;a=-0.4;
t=0:0.01:10;
ft=A*exp(a*t);
plot(t,ft);grid on;

%正弦信号
A=1;w0=2*pi;phi=pi/6;
t=0:0.01:8;
ft=A*sin(w0*t+phi);
plot(t,ft);grid on;

%抽样信号Sa(t)
t=-3*pi:pi/100:3*pi;
ft=sinc(t/pi);
plot(t,ft);grid on;
%or
sym t
Sa=sym(‘sin(t)/t’)
ezplot(Sa,[-10 10]) 
%三角波
y=sawtooth(2*pi*50*x,0.5);
hp=plot(x,y)

%锯齿波
y=sawtooth(2*pi*50*x,1);
hp=plot(x,y)

%方波
t = -1:0.001:1;
y = square(2 * pi * t,50);%周期是1,占空比是50%的方波
hp=plot(x,y)

%周期方波信号
plot(t,y), grid on;
%非周期三角波
y=tripuls(x,0.1);
hp=plot(x,y)
%非周期方波(也叫矩形信号）
t=0:0.001:4;
T=1;
ft=rectpuls(t-2*T,2*T);
plot(t,ft);grid on;axis([0 4 -0.5 1.5]);