昆明理工大学信号与系统B实验

灼华夫诸

已于 2024-06-13 11:17:23 修改

阅读量882

点赞数 34

文章标签： matlab 信号处理

于 2024-06-12 21:25:01 首次发布

本文链接：https://blog.csdn.net/2203_75285955/article/details/139634307

版权

前言

本文仅限昆明理工大学开设的信号与系统B2019版实验的MATLAB结果分享。文章部分内容可能有错误，仅代表个人观点。

实验一连续信号的时域描述与运算实验

1.上机实验前, 认真阅读实验原理, 掌握信号表示和信号运算的方法。

这部分将不过多讲述，认真阅读指导书即可。

2.利用 MATLAB 命令画出下列连续信号的波形图。

目标一： $(2-e^{-t})\varepsilon (t)$

目标二： $\varepsilon (cost)$

目标三：产生幅度为 1,周期为 1,占空比为 0.5 的周期矩形信号

由于这三个目标(子任务)之间没有关联性，因此不使用syms定义变量和函数(连续时间信号最好使用定义函数的方式)。此处利用离散时间信号的采样对函数进行赋值。

Fs = 100;    % 设置采样频率
T = 1 / Fs;  % 设置采样周期
t = -10:T:10;% 设置长度

% 目标2.1
f2_1 = (2-exp(-t)).*heaviside(t);

% 目标2.2
f2_2_cos = cos(pi*t/5);
f2_2 = heaviside(f2_2_cos);

% 目标2.3
for i = 1:length(t)
    if mod(i, Fs)  == 0
        f2_3(i) = 1;
    elseif mod(i, Fs * 2) == 1
        f2_3(i) = 0;
    else f2_3(i) = f2_3(i - 1);
    end
end

通过采样设置连续时间信号的弊端是：在采样频率较小采样周期过大时，输出结果(包括图像)的精度误差过大。最终得到的图像如下。

3.已知f(t)的波形如图所示,作出 f(t)+f(t)、f(t)·f(t)、f(t)的微分、f(t)的积分、f(3-4t),f(1-t/1.5)以及 f(t)*f(3-4t)的波形,并作出f(t)的奇、偶分量。

这里我将利用前文提到的，利用syms定义变量和函数，因为f(t)在这个任务中被多次利用与修改。从实验图可以看到，f(t)是一个分段函数。要设置分段函数，可以使用阶跃函数的加减，实现分段存在的0或1信号，在对应分段中.*(标量乘)函数就行了。

syms x;
f3 = @(x) (heaviside(x)-heaviside(x-1)).*x-heaviside(x-1)+heaviside(x-2);    % 函数的设置

% 目标3.1
f3_1 =@(x) f3(x) + f3(x);

% 目标3.2
f3_2 =@(x) f3(x) .* f3(x);

% 目标3.3
f3_3 = diff(f3, x); 

% 目标3.4
f3_4 = int(f3, x, 0, 2); 

% 目标3.5
f3_5_1 =@(x) f3(3 - 4 .* x);
f3_5_2 =@(x) f3(1 - x / 1.5);

% 目标3.6
f3_6 =@(x) conv(f3(x), f3(3 - 4 .* x));

% 目标3.7
f3_7_odd =@(x) (f3(x) - f3(-x)) / 2;
f3_7_eve =@(x) (f3(x) + f3(-x)) / 2;

f3_1对应两函数相加(图中目标一)，f3_2对应两函数的点乘(标量乘)(图中目标二)，f3_3对应f(t)的微分(图中目标三)，f3_4对应函数f(t)在[0,2]上的定积分(图中目标四)，f3_5_1和f3_5_2对应函数的尺度变换：f(3-4t)和f(1-t/1.5)(图中目标五和目标六)，f3_6对应两函数的卷积积分：f(t)*f(3-4t)(图中目标七)，f3_7_odd和f3_7_eve对应f(t)的奇、偶分量(图中目标八和目标九)。

实验二离散信号的时域描述与运算实验

1.上机实验前,认真阅读实验原理,掌握序列表示和运算的方法。

同上一个实验，这一部分不多描述，仔细阅读指导书即可。

2.利用 MATLAB 命令画出下列序列的波形图。

目标一： $(2-0.5^{-n})\varepsilon (t)$

目标二： $(\frac{3}{2})^{n}sin(\frac{n\pi }{5})$

离散时间信号的单位阶跃函数不建议使用heaviside，因为在零点的heaviside(0)是等于0.5的，在理论上连续时间信号的单位阶跃函数的零点确实是0.5，但是离散时间信号中的单位阶跃函数的零点应该是等于1的。为了代替函数heaviside，可以使用(x >= 0)，相当于是一个逻辑式，逻辑真为1，逻辑假为0。

离散时间信号也可以通过syms定义函数，在通过对函数进行采样得到。由于此处两个目标较简单，不使用对函数采样的方式描述。

x = 0:0.1:10;    % 设置区间
f2_1 =(2 - (0.5 .^ x)) .* (x>=0);
f2_2 =(1.5 .^ x) .* sin(pi .* x / 5);

离散时间信号的时间域一定是离散的，因此不能像函数的输出图像可以用平滑的曲线连接。

3.已知 f(n)的波形如图所示,作出f(n)的一阶后向差分、f(n)的累加和、 f(3-2n)以及 f(n)*f(3-2n)的波形。

这里我也没有利用对函数进行采样的方式，仍然是直接采样得到序列(涉及到序列变换的最好通过对函数的采样进行)。另外，序列一般不用f(t)表示，而利用x[n]表示。序列的差分类似于函数的微分，序列的累加类似与函数的积分。

% 设置必要参数
l = - 2;    % 起始
r = 5;      % 终止
l0 = 0;     % 有值起始
r0 = 4;     % 有值终止
s = l:1:r;  % 序列区间
s1 = l+1:1:r;% 差分结果区间
s2_1 = l/2 :1: r/2;
s2 = 3 - s2_1;% 变换后的区间

f3 = (s>=l0)-(s>=r0);

% 目标3.1
for i = 2:length(f3)
    f3_1(i-1) = f3(i) - f3(i-1);
end

% 目标3.2
f3_2(1) = f3(1);
for i = 2:length(f3)
    f3_2(i) = f3(i) + f3_2(i-1);
end

% 目标3.3
f3_3 =abs( (s2>= (3 - l0 / 2)) - (s2>= (3 - r0 / 2)) );

% 目标3.4
f3_4 = conv(f3, f3_3);
s4 = l * (l/2) :1: r * (r/2);

f3对应原序列f(n)，f3_1对应序列的一阶后向差分，f3_2对应序列的累加和，f3_3对应序列的变换f(3-2n)，f3_4对应序列的卷积和f(n)*f(3-2n)。由于网格的存在，因此我将线条的宽度加粗了。

实验三连续信号的频域分析实验

1.求图所示周期信号(T=2,T=1)的傅里叶级数,用 MATLAB 作出其前 3、21、45 项谐波的合成波形并与原信号作比较,作出其单边幅度谱和相位谱。

从实验三开始没有任务要求阅读指导书，但是最好也认真读一遍指导书。题目给的图是T=2时的波形，还要自己设置T=1时的波形，然后再进行傅里叶级数的展开、合成谐波、傅里叶变换。

Fs = 48; T = 1/Fs;    % 设置采样频率和采样周期
x = -6:T:6;l = length(x);    % 设置函数长度
for i = 1:(Fs / 2 + 1)
    fx2(i) = (i-1)/(Fs / 2);
    fx2(Fs + 2 - i) = (i-1)/(Fs / 2);
end    % 生成(T = 2)一个周期内的函数
for i = 0:1:5
    ft2((1 : (Fs / 2 + 1)) + i * 2 * Fs) = fx2((Fs / 2 + 1):(Fs + 1));
    ft2(((3 * Fs / 2 + 1):(Fs * 2 + 1)) + i * 2 * Fs) = fx2(1:(Fs / 2 + 1));
end    % 生成(T = 2)的波形
for i = 1:(Fs / 4 + 1)
    fx1(i) = (i-1)/(Fs / 4);
    fx1((Fs / 2 + 2) - i) = (i-1)/(Fs / 4);
end    % 生成(T = 1)一个周期内的函数
for i = 0:1:11
    ft1((1 : (Fs / 4 + 1)) + i * Fs) = fx1((Fs / 4 + 1):(Fs / 2 + 1));
    ft1(((3 * Fs / 4 + 1):(Fs + 1)) + i * Fs) = fx1(1:(Fs / 4 + 1));
end    % 生成(T = 2)的波形

n_3 = 3; n_21 = 21; n_46 = 46; n_f = Fs;    % 设置谐波阶次

ans1 = fft(ft1,l);
a_1_3 = ans1(1:n_3+1);
a1_3 = zeros(size(ans1));
a1_3(1:n_3+1) = a_1_3;
f_1_3 = real(ifft(a1_3));  
a_1_21 = ans1(1:n_21+1); 
a1_21 = zeros(size(ans1));
a1_21(1:n_21+1) = a_1_21;
f_1_21 = real(ifft(a1_21)); 
a_1_46 = ans1(1:n_46+1); 
a1_46 = zeros(size(ans1));
a1_46(1:n_46+1) = a_1_46;  
f_1_46 = real(ifft(a1_46));    % 对(T = 1)的函数进行谐波合成

ans2 = fft(ft2,l);
a_2_3 = ans2(1:n_3+1);
a2_3 = zeros(size(ans2));
a2_3(1:n_3+1) = a_2_3;
f_2_3 = real(ifft(a2_3));  
a_2_21 = ans2(1:n_21+1); 
a2_21 = zeros(size(ans2));
a2_21(1:n_21+1) = a_2_21;
f_2_21 = real(ifft(a2_21)); 
a_2_46 = ans2(1:n_46+1); 
a2_46 = zeros(size(ans2));
a2_46(1:n_46+1) = a_2_46;  
f_2_46 = real(ifft(a2_46));    % 对(T = 2)的函数进行谐波合成

a_1_f = ans1(1:n_f+1);
a1_f = zeros(size(ans1));
a1_f(1:n_f+1) = a_1_f;
f_1f = real(ifft(a1_f));  
a_2_f = ans2(1:n_f+1);
a2_f = zeros(size(ans2));
a2_f(1:n_f+1) = a_2_f;
f_2f = real(ifft(a2_f));    % 设置一个全波合成作为参照

a = abs(fft(ft1))/l;    % 通过快速傅里叶变换得到幅度谱
ph = angle(fft(ft1))/l; % 通过快速傅里叶变换得到相位谱

由于没有合理使用矩阵与向量的概念，因此这里的代码较长。

2.求图所示的单个三角脉冲(T=1)的傅里叶变换,并作出其幅度谱和相位谱

这个目标没有要求谐波的合成，因此直接进行快速傅里叶变换。

Fs = 10;T = 1 / Fs;        % 设置采样频率和采样周期
x = -6:T:6;l = length(x);  % 设置函数长度
fx = zeros(1,l);
for i = 1:(1 + 0.5 * Fs)
    fx(i + 5.5 .* Fs) = (i-1)/(0.5 .* Fs);
    fx(6.5 .* Fs - i + 2) = (i-1)/(0.5 .* Fs);
end    % 设置函数

a = abs(fft(fx))/l;    % 快速傅里叶变换得到幅度谱
ph = angle(fft(fx))/l; % 快速傅里叶变换得到相位谱

3.求不同占空比下, 周期矩形脉冲的幅度谱和相位谱,例如 $\frac{\tau }{T} = \frac{1}{8} , \frac{1}{4}$

Fs = 200;  
T = 1/Fs;
t = (0:Fs-1)*T;    % 函数的基本信息

dc_4 = 0.25; 
f_4 = zeros(size(t));  
for i = 1:length(t)  
    if (mod(t(i), FT) < dc_4 * FT) f_4(i) = 1;  
    else f_4(i) = 0; 
    end  
end    % 设置占空比0.25的波形

dc_8 = 0.125; 
f_8 = zeros(size(t));  
for i = 1:length(t)  
    if (mod(t(i), FT) < dc_8 * FT) f_8(i) = 1;  
    else f_8(i) = 0; 
    end  
end    % 设置占空比0.125的波形
a_4 = abs(fft(f_4))/Fs;    % 占空比0.25的幅度谱
ph_4 = angle(fft(f_4))/Fs; % 占空比0.25的相位谱
a_8 = abs(fft(f_8))/Fs;    % 占空比0.125的幅度谱
ph_8 = angle(fft(f_8))/Fs; % 占空比0.125的相位谱

类型相似的实验，不过多描述。由于是对比实验，因此将不同占空比的信号结果放在一起。

4.验证傅里叶变换的性质，如

目标一：时移性质:选取f(t)和 f(t-b), 幅频曲线相同, 只有相位不同

目标二：频移性质:选取 $\varepsilon (t)$ 和 $cos\Omega _{0}t\cdot \varepsilon (t)$ 或 $sin\Omega _{0}t\cdot \varepsilon (t)$

目标三：对称性质:选取 $Sa(\Omega _{0}t)$ 和 $g_{\tau }(t)$

目标四：尺度变换性: 选取 f(t)和 f(at)

目标四贴的代码是错误的，正确的代码私信我者发。不明白为什么任务的目标一有后半句的图像特征，而其他三个目标没有。一样的设置原函数和快速傅里叶变换，主要是对幅度谱和相位谱的对比分析。

% 时移
Fs = 200; b = 25; l = 1:1:Fs / 2 + 1;
oumiga = 1/4 * pi;
for i = l
    fo(i) = abs(Fs / 4 - abs(Fs / 4 + 1 - i));
    fb(i) = abs(Fs / 4 - abs(Fs / 2 - b + 1 - i));
end
ao = abs(fft(fo))/Fs;
ab = abs(fft(fb))/Fs;
pho = angle(fft(fo))/Fs;
phb = angle(fft(fb))/Fs;

% 频移
x = -Fs / 2 :1: Fs / 2;
for i = 1 : length(x)
    fe(i) = heaviside(x(i));
    fw(i) = cos(oumiga .* x(i)) .* heaviside(x(i));
end

ae = abs(fft(fe))/Fs;
aw = abs(fft(fw))/Fs;
phe = angle(fft(fe))/Fs;
phw = angle(fft(fw))/Fs;

% 对称
fsa = sin(x) ./ x;
for i = 1 : length(x)
    if (abs(x(i)) <= Fs / 4)
        fg(i) = 1;
    else fg(i) = 0;
    end
end
asa = abs(fft(fsa))/Fs;
ag = abs(fft(fg))/Fs;
phsa = angle(fft(fsa))/Fs;
phg = angle(fft(fg))/Fs;

% 尺度变换
k = 0.25;
fk = k .* fo;
ak = abs(fft(fk))/Fs;
phk = angle(fft(fk))/Fs;

这段程序不添加注释了，一样的注释内容。接下来的四张图依次是时移、频移、对称、尺度变换。

实验四离散信号的频域分析实验

1.对连续信号 $x_{a}(t)=Ae^{-jat}sin(\Omega _{0}t)\cdot \varepsilon (t)$ $(A=444.128, a=50\sqrt{2}\pi , \Omega _{0}=50\sqrt{2}\pi )$ 进行理想采样，可得到采样序列 $x(n)=x_{a}(nT)=Ae^{-jant}sin(\Omega _{0}nt)\varepsilon (t) 0<n<50$ 。图给出了 $x_{a}(t)$ 的幅频特性曲线，由图可以确定对 $x_{a}(t)$ 采用的采样频率。分别取采样频率为1kHz，300Hz和200Hz，画出所得采样序列x(n)的幅频特性 $\left | X(e^{j\omega }) \right |$ ，并观察是否存在频谱混叠。

老实说我并没有读懂这个实验的文本内容，依个人看法，就是用不同的采样频率得到不同的幅度谱，然后对三个幅度谱进行对比分析，没有相位谱的对比分析。

由于是对函数进行采样，可以先用syms定义函数，再对函数进行特定采样频率的采样(一般来说老师就希望这么做)。但是此处我仍然直接采样得到序列。

A = 444.128;  
a = 50 * sqrt(2) * pi;  
oumiga = 50 * sqrt(2) * pi; 
t = 0:0.0001:50;    % 设置函数特征
  
Xa = A * exp(-a * t .* 1i) .* (exp(1i * oumiga * t) - exp(-1i * oumiga * t)) ./ (2i);
Xa(t < 0) = 0;      % 设置函数

% 通过不同的采样频率得到不同的序列
Fs1 = 1000; Fs2 = 300; Fs3 = 200;  
n1 = 0:1/Fs1:50-1/Fs1; n2 = 0:1/Fs2:50-1/Fs2; n3 = 0:1/Fs3:50-1/Fs3;  
x1 = A * exp(-a * n1 .* 1i) .* (exp(1i * oumiga * n1) - exp(-1i * oumiga * n1)) ./ (2i);  
x2 = A * exp(-a * n2 .* 1i) .* (exp(1i * oumiga * n2) - exp(-1i * oumiga * n2)) ./ (2i);  
x3 = A * exp(-a * n3 .* 1i) .* (exp(1i * oumiga * n3) - exp(-1i * oumiga * n3)) ./ (2i);  

% 快速傅里叶变换得到幅度谱
aa = abs(fft(Xa))/length(t);
a1 = abs(fft(x1))/length(n1);
a2 = abs(fft(x2))/length(n2);
a3 = abs(fft(x3))/length(n3);

由于对题目理解的问题，得出的幅频谱与题目给出的幅频谱不一致。

2.设 x(n)=cos(0.48πn)+cos(0.52πn)

目标一：x(n)(0<n<10)时,求 x(n)的FFT 变换 X(k),并绘出其幅度曲线。

目标二：将 (1) 中的 x(n)以补零方式加长到 0≤n≤20, 求 X(k)并绘出其幅度曲线。

目标三：x(n)(0<n<100)时,求 X(k)并绘出其幅度曲线

目标四：观察上述三种情况下,x(n)的幅度曲线是否一致?为什么?

这个任务是变化序列的长度，没有改变采样频率和采样周期，并对比分析三个相同表达式不同长度序列的幅度谱。

这次我使用syms定义函数并取样，可以借鉴使用方法至之前的实验中。

syms n;
x = @(n) cos(0.48 .* pi .* n) + cos(0.52 .* pi .* n);    % 设置函数

n1 = 0:0.01:10; n2 = 0:0.01:20; n3 = 0:0.01:100;    % 设置区间
x1 = x(n1);  x2 = x(n2); x3 = x(n3);                % 得到序列

% 快速傅里叶变换得到幅度谱
a1 = abs(fft(x1))/length(n1);
a2 = abs(fft(x2))/length(n2);
a3 = abs(fft(x3))/length(n3);

延长了序列的定义域之后，对比发现幅频谱的高频区变化了。

3.（本任务没有主标题）

目标一：编制信号产生子程序,产生以下典型信号供谱分析用
$x1(n)=\left\{\begin{matrix} n+1 & 0\leq n\leq 3\\ 8-n & 4< n< 7\\ 0 & 其它n \end{matrix}\right.$ $x2(n)=cos(n\pi /4)$

$x3(n)=sin(n\pi /8)$

$x4(n)=cos8\pi t+cos16\pi t+cos20\pi t$

目标二：对信号 x1(n),x2(n),x3(n)进行两次谱分析,FFT 的变换区间 N 分别取 8 和 16，观察两次的
结果是否一致? 为什么?

目标三：连续信号 x4(n)的采样频率 fs=64Hz,N=16,32,64.观察 3 次变换的结果是否一致?为什么?

这个任务只说谱分析，所以应该放出幅度谱和相位谱，并对比分析。“变换区间N分别取8和16”这个也没有怎么理解，根据个人理解和目标三的定义，应该是说采样频率。

% 设置四个函数
syms n t;
x1_n =@(n) (n + 1) .* (heaviside(n) - heaviside(n - 3)) + (8 - n) .* (heaviside(n - 4) - heaviside(n - 7));
x2_n =@(n) cos(n .* pi ./ 4);
x3_n =@(n) sin(n .* pi ./ 8);
x4_t =@(t) cos(8 .* pi .* t) + cos(16 .* pi .* t) + cos(20 .* pi .* t); 

% 根据采样数进行序列的取样
n1 = 1:0.125:8; n2 = 1:0.0625:8;
x1_1 = x1_n(n1); x1_2 = x1_n(n2);
x2_1 = x2_n(n1); x2_2 = x2_n(n2);
x3_1 = x3_n(n1); x3_2 = x3_n(n2);

% 获取三个序列的幅度谱和相位谱
a11 = abs(fft(x1_1)); ph11 = angle(fft(x1_1));
a21 = abs(fft(x2_1)); ph21 = angle(fft(x2_1));
a31 = abs(fft(x3_1)); ph31 = angle(fft(x3_1));
a12 = abs(fft(x1_2)); ph12 = angle(fft(x1_2));
a22 = abs(fft(x2_2)); ph22 = angle(fft(x2_2));
a32 = abs(fft(x3_2)); ph32 = angle(fft(x3_2));

% 对第四个序列进行采样和获取谱分析
Fs = 64; T = 1 / Fs;
n1 = 16; n2 = 32; n3 = 64;
t1 = 0:(1 / n1):(n1 .* T);
t2 = 0:(1 / n2):(n2 .* T);
t3 = 0:(1 / n3):(n3 .* T);
x4_1 = x4_t(t1);
x4_2 = x4_t(t2);
x4_3 = x4_t(t3);
a41 = abs(fft(x4_1)); ph41 = angle(fft(x4_1));
a42 = abs(fft(x4_2)); ph42 = angle(fft(x4_2));
a43 = abs(fft(x4_3)); ph43 = angle(fft(x4_3));

分段函数的设置之前说过，利用单位阶跃函数的时移和相减，得到一个区间的1，类似门函数的时移。由此可以获取多个高1的门区间，再将这几个1的区间.*对应的函数，得到最后的分段函数。以下结果依次为：x1(n)、x2(n)、x3(n)、x4(n)四个原信号，前三个信号的谱分析，x4的谱分析。