使用matlab对方波信号进行仿真实验

实验任务

本实验来自于信号与系统五一大作业，实验任务如下：现有一输入信号为x_(t)为一方波信号，其波形如图所示。其中 T = 2 , T 1 = 0.5 T=2,T_1=0.5 T=2,T1​=0.5。用matlab对该信号进行如下仿真任务。

1. 用square()函数绘出该方波信号

2. 用谐波模拟，其中k取5，11，40

3. 绘制数a_k的图像，其中a_k为其傅里叶系数

4. 将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​∫T​x(t)​e−jkω0​tdx 。我们可以得到该方波信号的傅里叶系数展开式为

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ω0​T)​，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时a_k=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('合成的方波');

绘制数a_k的图像，其中a_k为其傅里叶系数

代码

%绘出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=+∞​ak​H（jkω0​）​ejkω0​t​

代码

%将用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绘制出的图像一致。