使用matlab对方波信号进行仿真实验
作者：互联网
使用matlab对方波信号进行仿真实验
实验任务
本实验来自于信号与系统五一大作业，实验任务如下：现有一输入信号为x(t)为一方波信号，其波形如图所示。其中
T
=
2
,
T
1
=
0.5
T=2,T_1=0.5
T=2,T1=0.5。用matlab对该信号进行如下仿真任务。
-
用square()函数绘出该方波信号
-
用谐波模拟，其中k取5，11，40
-
绘制数ak的图像，其中ak为其傅里叶系数
-
将x(t)输入到h（t）=e-tu(t)中，得到的信号为y(t) ,计算出y(t)的傅里叶系数并绘制其图像
用square()函数绘出该方波信号
代码
%用square()函数 仿真方波信号
t=-5:0.01:5;
x=0.5*(1+square(pi*(t+0.5)));
plot(t,x);
set(gca,'YLim',[0,1.5]);%改变y轴的坐标显示范围
**注释：**查阅官方文档我们可知，square()函数生成的方波信号的默认周期为2pi，范围为[-1,1],且在第一个周期即[0,2*pi]内在前半个周期为高，后半个周期为低，如下图所示。所以对函数进行如上变换。
图像
用谐波模拟
代码
%用5谐波模拟方波信号
t=-5:0.01:5;
T=2;
T1=0.5;
a0=2*T1/T;%直流部分
w=2*pi/T;
x=a0;
k0=5;%通过改变k0来改变谐波次数
num=1;
figure(1);
for k=1:1:k0
m=sin(k*w*T1)/(k*pi);
xk=2*m.*cos(k*w*t);%该式来自于奥本海姆《信号与系统》e.q3.31
x=x+xk;
subplot(k0,1,k);
xk=real(xk);%取其实部（实际上虚部为零，只是为了防止程序警告）
plot(t,xk);
title(['方波',num2str(k),'次谐波']);
end
figure(2);
x=real(x);%取其实部（实际上虚部为零，只是为了防止程序警告）
plot(t,x);
title('合成的方波');
图像
此时我们发现偶次谐波的数量级相对于奇次谐波的数量级很小，这是因为实际上奇次谐波是不存在的。根据傅里叶系数的展开公式 a k = 1 T ∫ T x ( t ) e − j k ω 0 t d x a_k=\frac{1}{T}\int_T\mathrm{x_{(t)}e^{-jk\omega_0t}\mathrm{d}x} ak=T1∫Tx(t)e−jkω0tdx 。我们可以得到该方波信号的傅里叶系数展开式为
a
k
=
{
s
i
n
(
k
ω
0
T
)
k
π
，
k
≠
0
2
T
1
T
，
k
=
0
a_k=\begin{cases} \frac{sin(k\omega_0T)}{k\pi}，k\neq0 \\ \frac{2T_1}{T}， k=0\end{cases}
ak={kπsin(kω0T)，k=0T2T1，k=0
将
ω
0
=
2
π
T
,
T
=
2
,
T
1
=
0.5
\omega_0=\frac{2\pi}{T},T=2,T_1=0.5
ω0=T2π,T=2,T1=0.5带入得
a k = { s i n ( π k / 2 ) k π , k ≠ 0 1 2 , k = 0 a_k=\begin{cases} \frac{sin(\pi k/2)}{k\pi},k\neq0 \\ \frac{1}{2}, k=0 \end{cases} ak={kπsin(πk/2),k=021,k=0
根据该式我们可以得知当k取偶数且k不为0时ak=0，所以我们得到偶次谐波为0。在matlab中由于自身计算算法的问题，所以是一个很小的数但不是0。
改进
我们对上述代码进行改进
%用5谐波模拟方波信号
t=-5:0.01:5;
T=2;
T1=0.5;
a0=2*T1/T;%直流部分
w=2*pi/T;
x=a0;
k0=5;%通过改变k0来改变谐波次数
if mod(k0,2)==1 %用于判断奇次谐波的个数
num=(k0+1)/2;
else
num=k0/2;
end
figure(1);
for k=1:1:num %只绘出奇次谐波
k1=2*k-1;
m=sin(k1*w*T1)/(k1*pi);
xk=2*m.*cos(k1*w*t);
x=x+xk;
subplot(num,1,k);
xk=real(xk);%取其实部（实际上虚部为零，只是为了防止程序警告）
plot(t,xk);
title(['方波',num2str(k1),'次谐波']);
end
figure(2);
x=real(x);%取其实部（实际上虚部为零，只是为了防止程序警告）
plot(t,x);
title('合成的方波');
绘制数ak的图像，其中ak为其傅里叶系数
代码
%绘出ak的图像
k0=40;
k1=-k0:1:k0;
ak=(sin(k1*pi/2))./(k1*pi);
[r,c]=size(ak);%读取行数和列数
for i=1:c
if abs(ak(i))<1e-10 %过于小的数实际上为0，这个我们在前面已经证明
ak(1,i)=0;
end
end
k2=0;
a0=0.5;%由于在原点处会出现无意义的值，不会绘制出来，所以在这里给补上
A=[ak,a0];
K=[k1,k2];
stem(K,A);
图像
将x(t)输入到h（t）=e-tu(t)中，得到的信号为y(t) ,计算出y(t)的傅里叶系数并绘制其图像
计算过程
首先计算H(S)
H
(
s
)
=
∫
−
∞
+
∞
h
(
τ
)
e
−
s
d
τ
=
∫
−
∞
+
∞
e
−
τ
u
(
τ
)
e
−
s
τ
d
τ
=
1
s
+
1
H_{(s)}=\int_{-\infty}^{+\infty}\mathrm{h_{(\tau)}e^{-s} d\tau} =\int_{-\infty}^{+\infty}\mathrm{e^{-\tau}u_{(\tau)}e^{-s\tau}d\tau} =\frac{1}{s+1}
H(s)=∫−∞+∞h(τ)e−sdτ=∫−∞+∞e−τu(τ)e−sτdτ=s+11
所以
y
(
t
)
=
∑
k
=
−
∞
k
=
+
∞
a
k
H
（
j
k
ω
0
）
e
j
k
ω
0
t
y_{(t)}=\sum_{k=-\infty}^{k=+\infty}a_kH_{（jk\omega_0）}e^{jk\omega_0t}
y(t)=∑k=−∞k=+∞akH（jkω0）ejkω0t
代码
%将用x输入h，输出信号为y
t=2:0.01:20;
T=2;
T1=0.5;
a0=2*T1/T;
w=2*pi/T;
y=a0;
k0=40;
for k=[-k0:1:-1,1:1:k0]
if mod(abs(k),2)=1 %仅取奇次谐波
ak=sin(k*w*T1)/(k*pi);
hs=1./(1j*k*w+1);
yk=ak*exp(1j*k*w*t).*hs;
y=y+yk;
end
end
plot(t,real(y)); %虚部为一些很小的数，不绘制
title('y的实部');
图像
没有绘制出虚部的原因是虚部是不存在的，可以证明y(t)为实数。
物理意义的探索
我们发现，该单位脉冲响应为一阶低通滤波器的单位脉冲响应，利用multisim可以进行仿真。在这里仅做定性的仿真
首先我们画出如下电路图
仿真结果如图所示
图形的形状与matlab绘制出的图像一致。
标签：仿真,ak,0.5,T1,谐波,k0,matlab,信号,pi 来源： https://blog.csdn.net/m0_46495731/article/details/119256826