信号与系统第三次实验

实验内容

题目一

对周期为 4 的方波信号𝑓(𝑡)进行谐波分析，使用多张子图分别绘制 1 至 3 次谐波、1 至 9 次谐波、1 至 33 次谐波、1 至 99 次谐波叠加时的图像，并在每个子图中同时画出原始方波信号𝑓(𝑡)的图像。绘图区间至少包括两个完整的方波周期并有完整标注。可以 手算求解傅里叶级数系数，也可以通过代码仿真计算。𝑓(𝑡)波形如图 1 所示。

**解：**该方波信号周期为4，且为奇函数，故$a_n=0$，所以

$$f(t)=\sum_{n=1}^{\infty} b_n \sin (n \omega t)$$

计算得到

$$b_n=\left\{\begin{array}{cc} 0 & n=2,4,6, \cdots \\ \frac{2}{n \pi} & n=1,3,5, \cdots \end{array}\right.$$

该方波的信号的傅立叶级数为

$$f(t)=\frac{2}{\pi}\left\{\sin (\frac{\pi}{2}t)+\frac{1}{3} \sin (\frac{\pi}{2}3 t)+\cdots+\frac{1}{2 n-1} \sin [\frac{\pi}{2}(2 n-1) t]+\cdots\right\} \quad n=1,2, \cdots$$

因此, 只要由 $b_n$ 计算出 $f(t)$ 各次谐波的振幅, 再根据各次谐波的频率, 即可利用 MATLAB 绘出周期方波的各次谐波叠加后的波形。

代码

m=50;
t = -4:0.01:4; %时域波形的时间范围-2~2，采样间隔 0.01 
n = round(length(t)/4); %根据周期方波信号的周期，计算 1/2 周期的数据点数
f = [0.5*ones(n,1);-0.5*ones(n,1);0.5*ones(n,1);-0.5*ones(n+1,1)]; %构造周期方波信号
y = zeros(m+1,max(size(t)));
y(m+1,:) = f';
x = zeros(size(t)); 
kk = 1; 
for k=1:2:2*m-1 %循环显示谐波叠加图形
 x = x+sin(pi/2*k*t)/k; 
 y((k+1)/2,:) = 2/pi*x; %计算各次谐波叠加和
 if k==3||k==9||k==33||k==99
    subplot(2,2,kk);
    plot(t,y(m+1,:)); 
    hold on; 
    plot(t,y((k+1)/2,:)); %绘制谐波叠加信号
    hold off; 
    grid; 
    axis([-4 4 -1 1]); 
    title(strcat('第1至',num2str(k),'次谐波叠加')); 
    xlabel('t','Fontsize', 8); 
    kk =kk+1; 
 end
end 

仿真图

题目二

求周期为 4 的信号𝑔(𝑡)的傅里叶级数系数{$𝑎_𝑘$}，并绘制𝑔(𝑡)的幅度谱、相位谱。可以手算求解傅里叶级数系数，也可以通过代码仿真计算。𝑔(𝑡)的波形如图 2 所示。

**解：**由图可知信号的周期为 $\mathrm{T}_0=4$, 所以 $\omega_0=\frac{2 \pi}{4}=\frac{\pi}{2}, g(t)$ 在区间 $[0,4]$ 内的表达式为

$$g(t)= \begin{cases}4(t-\frac{1}{2}) & 0< t \leq 1 \\ -1-t & -3 \leq t <0\end{cases}$$

求 $g(t)$ 的傅里叶级数系数$C_n$

$$C_n=\frac{1}{T_0} \int_{-3}^{1} g(t) e^{-j n \omega_0 t} d t=\frac{1}{4} [\int_{0}^{1} 4( t-\frac{1}{2}) e^{-j n\omega_0 t} d t+ \int_{-3}^{0} (-1-t) e^{-j n \omega_0 t} d t]$$

代码

clear;
clc;
close all;
syms c
dt = 0.001;
t = -4:dt:4;
T = 4;
w0 = 2 * pi / T;
n = -10:10;
Cn = zeros(size(n));
for N = -10:10
    ak=0;
    ak=1/T*(int((-1-c) * exp(-1i * N * w0 * c),c,-3,0)+int((4 * c - 2) * exp(-1i * N * w0 * c) ,c,0,1));
    Cn(N+11)=ak;
end
subplot(2,1,1);stem(n,abs(Cn));ylabel('Cn 的幅度谱');grid;
subplot(2,1,2);stem(n,angle(Cn));ylabel('Cn 的相位谱');grid;
xlabel('\omega/\omega_0');

仿真图

题目三

利用连续时间傅里叶级数的时移性质和时间反转性质求信号 𝑔(1 − 𝑡)的傅里叶级数系数{$𝑏_𝑘$}，并绘制𝑔(1 − 𝑡)的幅度谱和相位谱；使用{$𝑏_𝑘$}生成并绘制信号𝑔(1 − 𝑡)，以验证{$𝑏_𝑘$}的正确性。{$𝑏_𝑘$} 与{$𝑎_𝑘$}的关系应在报告中给出。

由连续时间傅里叶级数的时移性质和时间反转性质，在题目二的基础上更改

$$b_k=e^{-jk\omega_0(-1)}a_{-k}$$

频谱-相位图代码

clear;
syms c
dt = 0.001;
t = -4:dt:4;
T = 4;
w0 = 2 * pi / T;
N=input('N=');

n = -N:N;
Cn = zeros(size(n));
for k = -N:N
    ak=1/T*(int((-1-c) * exp(-1i * k * w0 * c),c,-3,0)+int((4 * c - 2) * exp(-1i * k * w0 * c) ,c,0,1));
    Cn(k+N+1)=ak;
end
Dn = zeros(size(n));
for i = 1:2*N+1
    Dn(i) = Cn(2*(N+1)-i)*exp(-1i*n(i)*w0);
end
subplot(2,2,1);stem(n,abs(Dn));ylabel('Dn 的幅度谱');grid on;axis([-10,10,0,0.7]);
subplot(2,2,2);stem(n,angle(Dn));ylabel('Dn 的相位谱');grid on;axis([-10,10,-4,4]);

subplot(2,2,[3 ,4]);
g = zeros(size(t));
for i = 1:2*N+1
    g = g + Dn(i) * exp(1i*(i-N-1)*w0 * t);
end
plot(t, g);grid on;

仿真图-N=50次的拟合情况

题目四

验证傅里叶级数的相乘性质。使用傅里叶级数的相乘性质求信号 ℎ(𝑡) = 𝑓(𝑡)𝑔(𝑡)的傅里叶级数系数{$𝑐_𝑘$}，并绘制ℎ(𝑡)的幅度谱和相位谱；使用{$𝑐_𝑘$}生成并绘制信号ℎ(𝑡)，以验证{$𝑐_𝑘$}的正确性。{$𝑐_𝑘$} 的推导过程应在报告中给出。

解： $f(t)$ 和 $g(t)$ 为两个周期为 $T=4$ 的信号, 且$f(t) \stackrel{\mathcal{F S}}{\leftrightarrow} a_k $，$g(t) \stackrel{\mathcal{F S}}{\leftrightarrow} b_k $

则

$$h(t)=f(t) g(t) \stackrel{\mathcal{F S}}{\leftrightarrow} c_k=a_k * b_k=\sum_l a_l b_{k-l}$$

证明: $\quad f(t) g(t) \stackrel{\mathcal{F S}}{\leftrightarrow} c_k=\frac{1}{T} \int_T f(t) g(t) e^{-j k \omega_0 t} d t$

$$\begin{aligned} c_k &=\frac{1}{T} \int_T\left(\sum_m a_m e^{j m \omega_0 t}\right)\left(\sum_n b_n e^{j n \omega_0 t}\right) e^{-j k \omega_0 t} d t \\ &=\frac{1}{T} \int_T\left(\sum_m \sum_n a_m b_n e^{j(m+n) \omega_0 t}\right) e^{-j k \omega_0 t} d t \\ &=\frac{1}{T} \sum_m \sum_n\left(\int_T a_m b_n e^{j(m+n-k) \omega_0 t} d t\right)=\sum_m \sum_n\left\{\begin{array}{c} 0, m+n \neq k \\ a_m b_n, m+n=k \end{array}=\sum_l a_l b_{k-l}\right. \end{aligned}$$

代码

syms c
dt=0.01;
t=-3:dt:1;

%f g h
f = 0.5*square(0.5*pi* t);
% plot(t, f);
% hold on;
g=(-1-t).*(-t>0)+4*(t-0.5).*(t>0)+(-1).*(t==0);
% plot(t, g);
% hold on;
h=f.*g;
subplot(221)
plot(t, h);ylabel('h(t)');grid on;

%系数
T = 4;
w0 = 2 * pi / T;
N=input('N=');
n = -N:N;
An = zeros(size(n));
Bn = zeros(size(n));
for k = -N:N
    bk=1/T*(int((-1-c) * exp(-1i * k * w0 * c),c,-3,0)+int((4 * c - 2) * exp(-1i * k * w0 * c) ,c,0,1));
    Bn(k+N+1)=bk;

    ak=1/T*(int(0.5* exp(-1i * k * w0 * c),c,-3,-2)+int(-0.5* exp(-1i * k * w0 * c),c,-2,0)+int(0.5* exp(-1i * k * w0 * c),c,0,1));
    % if k==0
    %     ak=0
    % else
    %     ak=-1i/(k*pi);
    % end
    An(k+N+1)=ak;
end

m=-2*N:2*N;
Cm=conv(An,Bn);
subplot(222)
stem(m,abs(Cm));axis([-10,10,0,0.3]);ylabel('频谱图');grid on;
subplot(223)
stem(m,angle(Cm));axis([-10,10,-4,4]);ylabel('相位图');grid on;

h1 = zeros(size(t));
for i = 1:4*N+1
    h1 = h1 + Cm(i) * exp(1i*(i-2*N-1)*w0 * t);
end
subplot(224)
plot(t, h1);grid on;ylabel(strcat(num2str(N),'次谐波叠加的h(t)')); 

仿真图-N=50次的拟合情况

思考题

连续时间傅里叶级数和离散时间傅里叶级数有何不同点，造成不同的原因是什么？

离散时间周期信号的傅里叶级数是有限项级数，不存在收敛问题也没有吉柏斯现象，而在连续时间周期信号情况下是一个无穷级数。

原因在于任何离散时间周期序列$x[n]$完全是由有限个参数（N个）来表征的，就是在一个周期内的N个序列值，傅里叶分析公式只是把这N个参数变换为一组等效的N个傅里叶系数值。相比之下一个连续时间周期信号在单个周期内有连续取值问题，这就要求用无限多个傅里叶系数来表示它，需要考虑项数趋于无穷多的求极限问题，这自然就产生了收敛问题。

收获与感想

通过本次实验我掌握了利用 MATLAB 实现对周期信号的频谱分析、实现信号的幅度调制的方法，直观地感受到傅里叶级数是如何进行展开与拟合地，加深了对周期信号傅里叶级数的理解。