matlab通过傅里叶级数生成周期三角波

引入

我们知道,任何周期函数都可以用正弦函数和余弦函数构成的无穷级数来表示。那么周期三角波如何通过matlab这一工具使用傅里叶级数相加的方法来表示呢?
我所选择生成的三角波的函数表达式为:

f(t) = t - [t] , [t]为不大于t的最大整数

它的函数图像为:在这里插入图片描述

原理

1. 复数表示

根据三角波的函数表达式:

f(t) = t - [t] , [t]为不大于t的最大整数

我们使用周期函数求傅里叶级数的公式:
a k = 1 T ∫ 0 T f ( t ) e − j k ω 0 t d ⁡ t a_k=\frac1T\int_0^Tf\left(t\right)e^{-jk\omega_0t}\operatorname dt ak=T10Tf(t)ejkω0tdt
此时我们可以使用傅里叶级数来表示三角波函数:
f ( t ) ^ = ∑ k = − N N a k e j k w 0 t \widehat{f(t)}={\textstyle\sum_{k=-N}^N}a_ke^{jkw_0t} f(t) =k=NNakejkw0t
其中N在数学公式中理应为正无穷+∞,但是在代码实现过程中为一个我们可以调整的有限变量值,它表示傅里叶级数相加的位数。它越大,三角波函数的生成情况越好。我们也可以通过从小到大调整这一变量值,来观察三角波函数生产情况的由坏变好。

2. 三角函数

我们知道,一个周期为2L的周期函数f(t)可以表示为:
f ( t ) = a 0 2 + ∑ n = 1 ∞ a n cos ⁡ ( n π L t ) + ∑ n = 1 + ∞ b n sin ⁡ ( n π L t ) f(t)=\frac{a_0}2+{\textstyle\sum_{n=1}^\infty}a_n\cos(\frac{n\pi}Lt)+{\textstyle\sum_{n=1}^{+\infty}}b_n\sin(\frac{n\pi}Lt) f(t)=2a0+n=1ancos(Lnπt)+n=1+bnsin(Lnπt)
其中:
a n = 1 L ∫ − L L f ( t ) cos ⁡ n π L t d ⁡ t a_n=\frac1L\int_{-L}^Lf(t)\cos\frac{n\pi}Lt\operatorname dt an=L1LLf(t)cosLnπtdt
b n = 1 L ∫ − L L f ( t ) sin ⁡ n π L t d ⁡ t b_n=\frac1L\int_{-L}^Lf(t)\sin\frac{n\pi}Lt\operatorname dt bn=L1LLf(t)sinLnπtdt
处理上述系数的方式有两种:(1)通过matlab计算,这种方法可以避免我们自己人工运算积分式子,但是会大大增加程序的计算量,增加程序运行时间;(2)通过人工计算出积分式子,得到系数的表达式,省去matlab中积分的运算。

若采用第二种方法,需将我们这次需要实现的周期三角波函数代入到上述公式,计算出:
a n = { 2 ∫ 0 1 t cos ⁡ ( 2 n π t ) d t = 0 n ≠ 0 1 n = 0 a_n=\left\{\begin{array}{lc}2\int_0^1t\cos(2n\pi t)dt=0&n\neq0\\1&n=0\end{array}\right. an={201tcos(2nπt)dt=01n=0n=0
b n = 2 ∫ 0 1 t sin ⁡ ( 2 n π t ) d t = − 1 n π bn=2\int_0^1t\sin(2n\pi t)dt=-\frac1{n\pi} bn=201tsin(2nπt)dt=nπ1
所以三角波函数可以表示为:
f ( t ) ^ = 1 + ∑ n = 1 ∞ − 1 n π sin ⁡ ( 2 n π t ) \widehat{f(t)}=1+{\textstyle\sum_{n=1}^\infty}{\textstyle\frac{-1}{n\pi}}\sin(2n\pi t) f(t) =1+n=1nπ1sin(2nπt)

matlab代码

复数

复数积分在matlab中的代码还没得到很好的解决

三角函数

方法1(不推荐,运行时间很长)

运行时间很长,当N输入为1000时,matlab运行时间约为20分钟左右。

clear;
x = -3:0.01:3;
N = input('please enter the N')
y = g(N);  %N为相加的项数
plot(x,y);

 %定义累加函数,其变量N决定相加项数的多少
function y=g(N) 
x = -3:0.01:3;
y = zeros(1,length(x))
y = a(0)/2+y;
for k = 1:N
ai= a(k).*cos(2*k*pi.*x)+b(k).*sin(2*k*pi.*x); 
y = y + ai;  %各项累加
end
end

%定义求an的函数
function ak = a(k)  
syms t;
ak = 2*int(t*cos(2*k*pi*t),t,0,1);
end

%定义求bn的函数
function bk = b(k)
syms t;
bk = 2*int(t*sin(2*k*pi*t),t,0,1);
end

方法2

直接代入an,bn的值,运行速度很快,当n输入为1000时,运行时间也仅为1s以内

clear
n = input("please enter the n"); %取无穷级数的前K项叠加来近似f(t)
t = -2:0.01:2;
for k = 2:n
    f(k,:)=(-1)/pi/(k-1)*sin(2*pi*(k-1)*t);
end
f(1,:) = 1/2;  %加入a0的值
F = sum(f,1);
plot(t,F);

运行结果

因为三角函数的方法一和方法二是等价的,所以上述两个程序运行产生的效果是相同的。随着改变n的大小,结果如下图:
在这里插入图片描述
可以很明显地发现:随着相加项数的增加,三角波生成情况越来越好。

相关推荐
©️2020 CSDN 皮肤主题: 游动-白 设计师:白松林 返回首页