MATLAB 信号处理仿真入门实验

MATLAB 信号处理仿真入门实验

实验目的：
• 熟悉 Matlab 工具的基本用法
• 掌握 Matlab 代码编写方法
• 理解序列的离散时间傅里叶变换
• 理解 DFT 结果的频谱能量泄露
• 理解 DFT 和 DTFT 的对应关系
• 理解信号加窗的作用
实验内容：
• 任务1、单音正弦信号采样序列的时域绘图
• 任务2、单音正弦信号采样序列的DFT结果绘图
• 任务3、有限长矩形窗序列的DTFT
• 任务4、有限长正弦序列的DTFT
• 任务5、 DTFT 和 DFT 关系探寻
• 最终任务：双音信号频率分析
正弦信号及其采样序列的数学描述如下：

仿真中需要设定的信号配置参数：
1. 信号采样率 fs
2. 信号幅度 Amp
3. 正弦信号频率 f
4. 正弦初始相位为 phy
5. 信号采样序列长度 N
任务1、单音正弦信号采样序列的时域绘图
• 按照以下两组配置参数修改代码
• 观察绘制出的信号时域图样
配置参数 1：
f1 = 1000; Amp1 = 1; phy1 = 0; f2 = 7000; Amp2 = 1; phy2 = 0;
配置参数 2：
f1 = 1000; Amp1 = 1; phy1 = 0; f2 = 9000; Amp2 = 1; phy2 = 0;
实验代码：

 % test with Matlab 2016
close all; 
clear;
fig_fname = 'sine_1_tone.jpg';
fs  = 8E3;  % 采样率
N   = 32;   % 向量长度
% 信号频率、幅度、初相位
f1  = 1000; Amp1 = 1; phy1 = 0;
f2  = 7000; Amp2 = 1; phy2 = 0;
% 信号向量的下标索引
n_idx = [0:N-1];
% 生成信号采样序列向量
x1 = Amp1*sin(2*pi*f1/fs*n_idx + phy1);
x2 = Amp2*sin(2*pi*f2/fs*n_idx + phy2);
% 生成绘图的标题文字
str_fs = num2str(fs); str_N = num2str(N);
str_f1   = num2str(f1); str_Amp1 = num2str(Amp1); str_phy1 = num2str(phy1);
str_f2   = num2str(f2); str_Amp2 = num2str(Amp2); str_phy2 = num2str(phy2);
% 合成绘图标题文字
title_str1 = ['f1:', str_f1,'Hz, ', 'fs:', str_fs,'Hz, ' 'N:', str_N, ', Amp1:' str_Amp1, ', Phy1:', str_phy1];
title_str2 = ['f2:', str_f2,'Hz, ', 'fs:', str_fs,'Hz, ' 'N:', str_N, ', Amp2:' str_Amp2, ', Phy2:', str_phy2];

% 绘图
h = figure;

subplot(2,1,1);
plot(n_idx, x1, 'color', 'blue'); % 连线方式绘图
grid on; hold; 
stem(n_idx, x1, 'color', 'red' ); % 离散样值方式绘图
title(title_str1, 'fontsize',14); % 绘制标题文字

subplot(2,1,2);
plot(n_idx, x2, 'color', 'blue'); % 连线方式绘图
grid on; hold;                                       
stem(n_idx, x2, 'color', 'red' ); % 离散样值方式绘图
title(title_str2, 'fontsize',14); % 绘制标题文字

仿真图：

配置参数2代码与配置参数1代码相似，只需修改具体数值，下面只给出绘图。

• 请回答: 数字信号处理的课本中说数字采样信号的归一化角频率 ω 的数值范围是宽度为 2π 区间，这是为什么？
归一化频率是信号频率除以采样频率的一半，如果将归一化频率转换为角频率，则将归一化频率乘以π ，所以归一化角频率是宽度为2π 的区间。
• 请自行上网搜索关键字“带通采样定理”思考其原理和可行性
带通采样定理：设带通信号m(t),其频率限制在fL与fH之间，带宽为B=fH-fL，如果最小抽样速率fs=2fH/m,m是一个不超过fH/B的最大整数，那么m(t),可以完全由其抽样值确定。
• 任务2、单音正弦信号采样序列的DFT结果绘图
• 离散傅里叶变换 DFT
– 拥有著名的快速算法 — FFT
– DFT的主要用途：分析信号频谱
• 本实验的内容 – 计算确知信号（正弦信号）采样序列的DFT谱
– 设定不同的仿真配置，观察正弦信号的DFT谱
• 改变DFT的计算长度
• 改变输入正弦信号采样序列的频率
• 完成实验后思考问题 – 序列的DFT结果和其连续版本信号的频谱完全一致么？
实验代码：
主函数

% file: sine_1_tone_dft.m
% test with Matlab 2016
close all;
clear;
fig_fname_0 = 'sine_1_tone_dft_0.jpg';
fig_fname_1 = 'sine_1_tone_dft_1.jpg';
fig_fname_2 = 'sine_1_tone_dft_2.jpg';
fs  = 8E3; 
f0_0  = 500; f0_1  = 600; f0_2  = 500; 
N_0   = 32 ; N_1   = 32 ; N_2   = 35 ;

% plot each case
func_sine_1_tone_dft_plot(f0_0, fs, N_0, fig_fname_0);
func_sine_1_tone_dft_plot(f0_1, fs, N_1, fig_fname_1);
func_sine_1_tone_dft_plot(f0_2, fs, N_2, fig_fname_2);

子函数

% test with Matlab 2016
% file: func_sine_1_tone_dft_plot.m
% 绘图程序：绘制信号采样序列的 时域波形，
% 线性尺度DFT幅度谱，对数尺度DFT幅度谱
function func_sine_1_tone_dft_plot(f0, fs, N, fig_fname);
  n_idx = [0:N-1];  % 信号序列的下标向量 
  x = sin(2*pi*f0/fs*n_idx);    % 生成正弦信号序列
  y = abs(fft(x)); y_dB = 20*log10(1E-6+y); % 计算 DFT 结果线性幅度谱及对数幅度谱

  h = figure;
  % 绘图横轴左右边界
  x_left = -1; x_right = N+1;
  % 时域绘图 
  subplot(3,1,1);
  stem(n_idx, x); grid on; hold;
  plot(n_idx, x); xlim([x_left, x_right]);
  xlim([x_left, x_right]);
  title_str = ['f0:', num2str(f0),'Hz, ', 'fs:', num2str(fs),'Hz, ' 'N:', num2str(N)];
  title(title_str, 'fontsize',14);
  % 绘图线性幅度谱
  subplot(3,1,2); stem(n_idx, y); grid on;
  xlim([x_left, x_right]);
  title_str = ['DFT Magtitude Linear Scale, '];
  title_str = [title_str, ' Max:', num2str(max(y)), ', Min:', num2str(min(y))];
  title(title_str, 'fontsize',14);
  % 绘图对数幅度谱
  subplot(3,1,3); stem(n_idx, y_dB);  grid on;
  xlim([x_left, x_right]);
  title_str = ['DFT Magtitude dB Scale, '];
  title_str = [title_str, ' Max:', num2str(max(y_dB)), ', Min:', num2str(min(y_dB))];
  title(title_str, 'fontsize',14);

  print(h, '-djpeg', fig_fname);  % 保存绘图结果文件
end

• 实验代码说明
– sine_1_tone_dft.m ： 实验入口代码
– 函数说明
• function func_sine_1_tone_dft_plot(f0, fs, N, fig_fname);
• 功能：绘制正弦信号采样序列的DFT幅度谱，然后保存曲线图
• 参数 f0：正弦信号频率
• 参数 fs ：采样率
• 参数 N ：序列长度
• 参数 fig_fname：保存的绘图文件名（jpg格式）
仿真图：

采样序列中包含了2个完整的信号周期的采样值，DFT幅度谱分析结果有2根峰值谱线，其余位置上均是0。

采样序列中包含了2个完整的信号周期，以及 一个不完整信号周期的采样值，DFT幅度谱分析结果有2根极大值谱线，但是其余位置上不为0。

采样序列中包含了2个完整的信号周期，以及 一个不完整信号周期的采样值，DFT幅度谱 分析结果有2根极大值谱线，但是其余位置上不为0。
总结规律
– 在采样序列包含整数个信号周期的情况下，信号采样后的DFT分析结果和连续信号的傅里叶变换频谱最相似。
–频谱泄露的最根本原因还在于信号的非周期截断，可以从时域和频域两方面来理解：
(1) 从时域上，傅里叶变换的潜在假设为待处理的有限信号为周期性无限信号的周期主体，即假设原始信号为当前有限信号的无限个周期延拓。因此，举个简单例子。我们截取50HZ 正弦信号的一个周期，其无限延拓就是最原始的50HZ正弦信号，因此一个周期的有限信号即可代表其原始的无限信号；若截取的有限信号不是50HZ信号的整数倍周期，可知该有限信号的无限延拓不可完全的复原原始的50HZ无限信号，其首尾连接处出现断续，从而引入高次谐波分量，产生频谱泄露。
(2)从频域上，假设信号周期T，频率F， 采样周期Ts，采样频率Fs， 采样点数N，则整数周期截断的物理表达为N*Ts = m*T <=> N*1/Fs=m*1/F <=> F=m*Fs/N. 在频域上，Fs对应2π，且频域分辨率为2π/N。 因此F=m*Fs/N意义为信号频率在FFT后的第m根谱线上。反之，当非周期截断时，则无法满足F=m*Fs/N，信号的频率成分分散k*2π/N的频率点上。
由于工程实际处理的信号基本为非平稳多谐波信号，因此频谱泄露不可避免。为了改善频谱泄露，可行的方向包含以下几点：
1. 增大FFT变换的点数N。 通过增大N，一方面提高频域分辨率，更大可能满足F=m*Fs/N， 另一方面，压缩sinc的瓣宽，降低泄露水平。
2. 选用合适的窗函数。根据不同的需求来选择不同特性的窗函数，主瓣宽但旁瓣衰减大的窗或是主瓣窄但旁瓣相对衰减小的窗。
• 任务3、有限长矩形窗序列的DTFT
背景知识：DTFT的重要性
• 计算机的局限性
– 只能计算有限长的数据（无限长的数据，只存在于理论分析的世界里）
– 只能计算离散化的数据（连续变量的函数，计算其函数值时，只能 在离散化的自变量的抽样点上进行）
• 离散时间序列傅里叶变换（DTFT）的重要性
– 在无穷长的尺度上累加
– 输出的结果函数是连续变量函数
– 确定的离散序列，无论有限长、无限长
• 其DTFT变换结果，具有理论唯一性
– 理论唯一性的重要在于：
• 可以使用多种工程计算方法无限逼近
• 可以用来检验工程计算方法的正确性和有效性
实验代码：
主函数：

% file: plot_dtft.m
% test with Matlab 2007
clc 
clear
close all

w_min = -1*pi;
w_max =  1*pi;

w_delta = pi/1000;

n_min =   0;
n_max =   3 ;


n = [n_min:1:n_max];
w = [w_min:w_delta:w_max];

N = length(n);

x = ones(1,N); 
x_dscp_str = 'x(n)';
Y = func_calc_dtft(x, n, w); 

h = func_plot_dtft(x,n,w,Y, x_dscp_str);

子函数1：

% file: func_plot_dtft.m
% test with Matlab 2007
% Y[w] = DTFT(x[n]);
% plot time sequence x[n], real/imag part of Y[w], abs of Y[w]
function h = func_plot_dtft(x,n,w,Y, x_dscp_str);
h = figure; 

N =