信号处理仿真基础实验

Matlab基本操作

在本部分内容，我们先熟悉一下Matlab的基本操作，这有助于后面的使用它进行信号的合成与分析。

生成等差序列行向量和列向量

生成全零和全一的向量和矩阵

向量元素的寻址

向量部分的寻址

查阅帮助文档

标量和向量的乘积

向量和向量的对应位置相乘

向量作为函数的参数

for 循环

if 分支判断

函数与作图

清除变量和关闭全部绘图窗口

练习作业 向量的寻址赋值

音频波形信号波形生成 及分析

生成正弦信号的采样序列


$$\begin{equation} y = A \cdot \sin(2\pi f_0 t + \phi) , ~~ t = \frac{n}{f_s}, n = 0 \cdots N-1 \label{eq:sine_t} \end{equation}$$

信号保存为WAV文件

N = 24000; fs = 8000; f0 = 400; phy = 0; A = 0.1;
filename = 'sine.wav';
t = [0:N-1]/fs;
y = A*sin(2*pi*f0*t+phy);
audiowrite(filename,y,fs);

相位与频率的关系初探：生成线性扫频信号的采样序列


$$\begin{equation} y = A \cdot \sin(\phi(t) + \phi_0) \\ f(t) = \frac{d\phi(t)}{dt} = \frac{f_1 - f_0}{(N-1)/f_s} \cdot t + f_0 \\ y = A \cdot \sin(2\pi \int_{0}^{t}f(\theta)\cdot d\theta + \phi_0) \\ y = A \cdot \sin(2\pi \frac{1}{2}\frac{f_1 - f_0}{(N-1)/f_s} \cdot t^2 + f_0\cdot t + \phi_0) \\ t = \frac{n}{f_s}, n = 0 \cdots N-1 \label{eq:chirp_t} \end{equation}$$

数字角频率的含义

背景介绍

实验 1


N = 100; fs = 100; f0 = 10; f1 = 110;  ;
t = [0:N-1]/fs;
x1 =  sin(2*pi*f0*t );
x2 =  sin(2*pi*f1*t );
figure;
subplot(2,1,1);hold;stem(x1, 'filled');title('x1(n)');plot(x1); xlabel('n');
subplot(2,1,2);hold;stem(x2, 'filled');title('x2(n)');plot(x2); xlabel('n');

实验 2


t_delta = 0.001;  %  time precision  of  analog wave plotting , in seconds

fs = 100;   % sample rate, in Hz
f1 = 10  ;   % signal 1 frequency , in Hz    
f2 = 110 ;   % signal 2 frequency , in Hz    

t_len  = 0.2  ; % time length of wave plotting

t_idx = [0:t_delta:t_len]; % time index

% get analog wave x1(t) and x2(t)
x_1_t = sin(2*pi*f1*t_idx); 
x_2_t = sin(2*pi*f2*t_idx); 

% get the sampled digital sampled serials x1(n) and x2(n)
Ts = 1/fs; % sample period
t_Ts_idx = [0:Ts: t_len];
x_1_n = sin(2*pi*f1*t_Ts_idx);
x_2_n = sin(2*pi*f2*t_Ts_idx);

% plot analog wave and digital serials in the same figure
close all
figure;
subplot(3,1,1);
h1 = plot(t_idx*1E3, x_1_t,'blue', t_idx*1E3, x_2_t, 'red');
legend(h1,'x1(t)','x2(t)','Location','northeast');hold; 
stem(t_Ts_idx*1E3, x_1_n,   'blue', 'Marker','O','MarkerSize',6,'LineWidth',1.5,'LineStyle',':');
stem(t_Ts_idx*1E3, x_2_n, 'red',  'Marker','X','MarkerSize',12,'LineWidth',1.5,'LineStyle','--');
title(['x1(t) ',int2str(f1),' Hz, ','x2(t) ',int2str(f2),' Hz,',' fs ',int2str(fs),' Hz']), xlabel('time(ms)');

% plot digital serials x1(n) and x2(n) in different figure
subplot(3,1,2);
stem(t_Ts_idx*1E3, x_1_n, 'filled', 'blue');title(['x1(n) ',int2str(f1),' Hz']);xlabel('time(ms)');
subplot(3,1,3);
stem(t_Ts_idx*1E3, x_2_n, 'filled', 'red ');title(['x2(n) ',int2str(f2),' Hz']);xlabel('time(ms)');

实验 3

信号的频域分析

有限长信号序列的离散傅里叶变换（DTFT） 的解析计算

本节，我们关注一个具有最简单频谱的信号，余弦信号 x(t) = cos(Ω0t)，经过采样后的有限长序列的频谱形式。


$$\begin{equation} X({j\Omega}~) = \int_{-\infty}^{\infty} x(t)~e~^{-j\Omega t} dt \label{eq:FT_def} \end{equation}$$

clear; close all; clc;
Ts              = 0.1   ;   
Omega_span      = 10*pi ;    
Omega_precison  = 0.01  ;   
N               = 36    ;
Omega_0         = 0.5*2*pi; 
Omega_vec       = [-Omega_span/2: Omega_precison: Omega_span/2];
asinc=@(Omega,N,Ts) sin(1/2*Omega*(2*N+1)*Ts)./sin(1/2*Omega*Ts); 
Omega_vec_n = Omega_vec - Omega_0; Omega_vec_p = Omega_vec + Omega_0;
X_rs_jOmega_n = pi*asinc(Omega_vec_n, N, Ts);
X_rs_jOmega_p = pi*asinc(Omega_vec_p, N, Ts); 
X_rs_jOmega   = X_rs_jOmega_n + X_rs_jOmega_p;
t_vec = [0:N-1]*Ts; 
x_rs  = cos(Omega_0*t_vec);
figure;
subplot(4,1,1); plot(t_vec, x_rs)  ;title('cos(\Omega t)');xlabel('time(sec)');grid on
subplot(4,1,2); plot(Omega_vec/pi, X_rs_jOmega   );title('X(j\Omega)  ' );xlabel('\Omega /\pi (rad/\pi)');grid on
subplot(4,1,3); plot(Omega_vec/pi,  X_rs_jOmega_n);title('X(j\Omega) N' );xlabel('\Omega /\pi (rad/\pi)');grid on
subplot(4,1,4); plot(Omega_vec/pi,  X_rs_jOmega_p);title('X(j\Omega) P' );xlabel('\Omega /\pi (rad/\pi)');grid on

练习作业 有限长矩形窗序列的频域特征

练习作业 有限长余弦序列的频域特征

有限长信号序列的离散傅里叶变换（DTFT） 的 编程计算

　　 离散序列的DTFT计算式，目前的主流课本中，大部分是直接给出定义式，为了加强离散序列的信号物理意义， 本节我们从一个连续信号的采样序列的傅里叶变换的角度，导出DTFT的定义式，从而强调该定义式中，各个变量在连续信号域的物理意义。

DTFT计算表达式的导出

上一部分内容中，我们讨论并用Matlab代码绘制了有限长余弦信号采样序列的频谱 Xrs(jΩ) 。 但是 Xrs(jΩ) 的计算过程是使用解析的方法人工推导出的公式 $\eqref{eq:seq_Xf}$的结果 。

　　解析推导的方法虽然精确， 但是仅适用于具有简单表达式的规则信号， 本节我们根据定义式， 采用计算机编程近似计算信号的有限长采样序列的频谱。

　　对于连续时间信号 xc(t) 经过 采样率 fs 即 采样周期 Ts = 1/fs 抽样后， 通过加矩形窗（rectangular window）选取从时刻 t = 0 到时刻 t = (N − 1)Ts ， 共 N 个样值构成抽样序列信号 xrs(t) ， 根据公式 $\eqref{eq:seq_xn}$有
$$ \begin {equation}\begin{aligned} \\ x_{rs}(t) &= x_c(t)\cdot\sum_{n = 0} ^ {N-1} \delta[t - n T_s] \\ &= \sum_{n = 0} ^ {N-1} x_c(nT_s) \cdot \delta[t - n T_s] \end{aligned} \label{eq:seq_xrst_N} \end{equation} $$

则信号xrs(t) 的频谱 Xrs(jΩ) 可按照以下公式计算
$$ \begin {equation}\begin{aligned} X_{rs}(j\Omega) &= \int_{-\infty} ^{+\infty} \sum_{n = 0} ^ {N-1} x_c(nT_s) \cdot \delta[t - n T_s] e^{-j\Omega t}dt \\ &= \sum_{n = 0} ^ {N-1} x_c(nT_s)\cdot e^{-j~\Omega~nT_s } \\ &= \sum_{n = 0} ^ {N-1} x_c(nT_s)\cdot e^{-j~2\pi~(\Omega/\Omega_s)~n }, ~~~\Omega_s = \frac{2\pi}{T_s} \end{aligned} \label{eq:Rs_Omega_nTs} \end{equation}$$

从公式 $\eqref{eq:Rs_Omega_nTs}$ 可见：

由此， 我们只需要研究xrs(t) 在一个周期内的特性， 就可以掌握其全部信息


$$ \begin {equation}\large\begin{aligned} X_{rs}(j\omega \cdot f_s) &= \sum_{n = 0} ^ {N-1} x_c(nT_s)\cdot e^{-j~\omega~n } \end{aligned} \label{eq:Xrs_j_omega_Ts} \end{equation}$$

记 Xrs(jω ⋅ fs) 为 Xrs(e jω) ， 以及 xc(nTs) 为 x[n] 则 式 $\eqref{eq:Xrs_j_omega_Ts}$ 变为：
$$ \begin {equation}\large\begin{aligned} X_{rs}(e^{~j\omega}) &= \sum_{n = 0} ^ {N-1} x[n]\cdot e^{-j~\omega~n } \end{aligned} \label{eq:X_e_j_omega_x_n} \end{equation}$$

观察式 $\eqref{eq:X_e_j_omega_x_n}$ ， 该式正是 《数字信号处理》 课本中定义的有限长序列的“离散时间序列傅里叶变换（DTFT）” ， 回顾式 $\eqref{eq:seq_xrst_N}$, $\eqref{eq:Rs_Omega_nTs}$, $\eqref{eq:Xrs_j_omega_Ts}$ ，其推导过程简介如下：

有限长信号序列的 DTFT 计算方法

式 $\eqref{eq:X_e_j_omega_x_n}$ 中，频域变量 ω 是一个计算机无法精确表示 的连续变量，但是我们可以计算出 ω 变量在 [0, 2π] 区间上的等 间隔采样位置处的频谱函数 X(ejω) 的数值。

于是，DTFT的计算过程可以用向量表示如下： 设向量$\overrightarrow{\omega} = [0,~\Delta\omega,~2\Delta\omega, \cdots, (K-1)\Delta\omega]^T$， $\overrightarrow{n} = [0,1,2,\cdots, N-1]$ ， $\overrightarrow{x} = [x(0),x(1),x(2),\cdots,x(N-1)]^T$，以及 频谱向量 $\overrightarrow{X}(e^{j\omega})$ ，则有：


$$\begin {equation}\large\begin{aligned} \overrightarrow{X}(e^{j\omega}) = e^{-j\cdot \overrightarrow{\omega} \cdot \overrightarrow{n}} \cdot \overrightarrow{x} \end{aligned} \label{eq:X_w_vector} \end{equation}$$

则


$$\begin{equation*} \overrightarrow{X}(e^{j\omega})= {\Huge{e}}^ {-j \begin{bmatrix} \Delta\omega \\ 2\Delta\omega \\ \vdots \\ (K-1)\Delta\omega \end{bmatrix} \cdot \begin{bmatrix} 0 & 1 & \cdots & N-1 \end{bmatrix} } \cdot \begin{bmatrix} x(0) \\ x(1) \\ \vdots \\ x(N-1) \end{bmatrix} \end{equation*}$$
即：
$$\begin{equation*} \overrightarrow{X}(e^{j\omega})= {\Huge{e}}^{-j \begin{bmatrix} 0 & \Delta\omega & 2\Delta\omega&\dots & (N-1)\Delta\omega \\ 0 & 2\Delta\omega & 4\Delta\omega&\dots & 2(N-1)\Delta\omega \\ 0 & 3\Delta\omega & 6\Delta\omega&\dots & 3(N-1)\Delta\omega \\ \vdots &\vdots &\ddots &\vdots &\vdots \\ 0 & (K-1)\Delta\omega & 2(K-1)\Delta\omega&\dots & (N-1)(K-1)\Delta\omega \\ \end{bmatrix} } \cdot \begin{bmatrix} x(0) \\ x(1)\\ \vdots \\ x(N-1) \end{bmatrix} \end{equation*}$$

close all; clear; 

dtft_calc=@(xn_col_vec,w_delta) ... 
    exp(-j*(([0:w_delta:2*pi].')*[0:1:length(xn_col_vec)-1])) * xn_col_vec;

N   = 100; fs = 100; f0 = 11; phy =7*pi/10; 
w_delta = 0.01*pi;
t   = ([0:N-1]).'/fs;
xn  = cos(2*pi*f0*t + phy);
% xn  = ones(N,1);

yn = dtft_calc(xn, w_delta);
w_idx = [0:w_delta:2*pi];

[max_abs,nn] = max(abs(yn(1:floor(length(yn)/2))));
max_w_pi = w_idx(nn)/pi;
max_ang = angle(yn(nn))/pi;

figure;
title_str =  ['xn: N   = ', num2str(N) ', fs = ', num2str(fs) , 'Hz, f0 = ', num2str(f0),'Hz,  phy =' , num2str(phy/pi), 'pi']; 
freq_est_str = [ ' w-delta = ', num2str(w_delta/pi),'pi,',  ' ESTIMATED: MAX Modular, omega is : ', num2str(max_w_pi),'pi', ', That is :',num2str(max_w_pi*fs/2), ' Hz', ', amp value is : ', num2str(max_abs), ', Angle is : ', num2str(max_ang),'pi']

subplot(5,1,1);stem([0:N-1], xn,'filled'     ); title(title_str); grid on;
subplot(5,1,2);    plot(w_idx/pi, abs(yn))  ; title(['Modular: ', freq_est_str] ); xlabel('\omega/ \pi : rad/ \pi');grid on;
subplot(5,1,3);    plot(w_idx/pi, angle(yn)/pi); title('Angle  ');xlabel('\omega/ \pi : rad/ \pi'); ylabel('\Phi/ \pi: rad/\pi'); grid on;
subplot(5,1,4);    plot(w_idx/pi, real(yn)) ; title('Real   '); xlabel('\omega/ \pi : rad/ \pi');grid on;
subplot(5,1,5);    plot(w_idx/pi, imag(yn)) ; title('Imag   '); xlabel('\omega/ \pi : rad/ \pi');grid on;

disp([title_str, freq_est_str]);

练习作业，正弦信号的 DTFT 结果分析

离散傅里叶变换（DFT）和离散时间序列傅里叶变换（DTFT）的关系

　　离散傅里叶变换（DFT） 是信号处理领域里最重要的正交变换， 其快速算法被称为 “快速傅里叶变换（FFT）”。

$\require{extpfeil}\Newextarrow{\rightleftharpoons}{5,5}{0x2194}$
$$ \begin {equation}\large\begin{aligned} \left\{\begin{array}{ll} x[n] &\rightleftharpoons[]{ \mathcal{DTFT}} & X(e^{j~\omega}) &: & X(e^{j~\omega}) &=&& \sum\limits_{n = 0}^{N-1}x[n]\cdot e^{-j~\omega~n } \\x[n] &\rightleftharpoons[]{ \mathcal{DFT}} & X[k] &: & X[k] &=&& \sum\limits_{n = 0}^{N-1}x[n]\cdot e^{-j~(2\pi/N)~nk }, \quad 0 \leq k \leq N-1 \end{array}\right. \end{aligned} \label{eq:DTFT_vs_DFT} \end{equation}$$

　　观察公式$\eqref{eq:DTFT_vs_DFT}$ 可以发现， DTFT和DFT 的区别有：

　　 经过观察可以发现 :

练习作业，正弦信号的 DFT 结果分析

练习作业，正弦信号的补零 DFT 结果分析

窗函数特性

　　计算机只能处理有限长的信号数据序列，有限长的信号序列等效于对无限长的序列在时域乘以有限长的矩形窗函数。 在之前的内容中，我们从模拟频率域的视角，分析了有限长矩形窗的采样信号的频谱 Rs(jΩ)， 然后又从连续信号的有限长抽样序列的频谱Xrs(jΩ)进行尺度变换得到了信号的归一化角频率频谱Xrs(jω)， 即：X(e jω) 　　本节我们将着重从数字角频率域研究问题，首先我们研究不同长度的矩形窗的频谱特性，然后再研究其他形状的窗函数的特性。

矩形窗函数

　　对于矩形窗而言，设其第一个分量值从 n = 0 位置开始， 则其唯一的参数为窗的长度 N 。则其 DTFT变换可以表示如下：


$$ \begin{equation}\large\begin{aligned} \\ r[n] &= \left\{\begin{matrix} 1,~~ 0 \leq n< N \\ 0,~~ n \geq N \end{matrix}\right. \\ R(e~^{j\omega~}) &= \mathcal{DTFT}\{r[n]\} = \sum\limits_{n = 0}^{N-1} e^{-j~\omega~n } \\ &= \frac{1 - e^{j\omega N}}{1 - e^{j\omega}} \\ &= \frac{e^{-j\omega N/2}[e^{-j\omega N/2} - e^{j\omega N/2}]}{e^{-j\omega /2}[e^{-j\omega/2} - e^{j\omega/2}]} \\ &= e^{-j\omega(N-1)/2}\cdot\frac{\sin(\omega N/2)}{\sin(\omega/2)} \label{eq:rect_DTFT} \end{aligned} \end{equation} $$

实验代码：矩形窗函数频域分析

close all; clear; clc;
dtft_calc=@(xn_col_vec,w_delta) ... 
    exp(-j*(([0:w_delta:2*pi].')*[0:1:length(xn_col_vec)-1])) * xn_col_vec;

N = 10;

xn = ones(N,1);
n_idx = [0:N-1].';
w_delta = 0.0001*pi;
yn = dtft_calc(xn, w_delta);
yn = fftshift(yn);
yn_dB = 20*log10(abs(yn));
yn_dB_norm  = yn_dB - max(yn_dB);
w_idx = [-pi:w_delta:pi]/pi;
figure;  
subplot(2,2,[1 2]);stem(n_idx,xn);hold;plot(n_idx,xn);
subplot(2,2,3);plot(w_idx,abs(yn)); grid on;
subplot(2,2,4);plot(w_idx,yn_dB_norm); grid on;ylim([-100,0]);

实验作业

使用以上代码，修改参数 N 运行代码回答问题

几种典型窗函数

　　历史上，人们在使用傅里叶变换分析统计数据的周期性规律时，发现先对数据进行加窗再进行傅里叶变换可以得到不同的分析效果。 有几种著名的窗函数被设计出来，分别具有不同的主瓣宽度和旁瓣抑制比。我们需要对这几种窗函数建立起定性的认识，因为 它们经常被用于各种信号分析工具当中。

汉宁窗 (Hanning window)


$$\begin {equation}\begin{aligned} \\w(n)= \left\{\begin{array}{ll} 0.5−0.5\cos(2\pi \frac{n}{N-1}), & 0\leq n < N \\ 0, & n \geq N \end{array}\right. \end{aligned} \label{eq:Hanning_win_def} \end{equation}$$

汉明窗 (Hamming window)


$$\begin {equation}\begin{aligned} \\w(n)= \left\{\begin{array}{ll} 0.54−0.46\cos(2\pi \frac{n}{N-1})), & 0\leq n < N \\ 0, & n \geq N \end{array}\right. \end{aligned} \label{eq:Hamming_win_def} \end{equation}$$

巴特莱特窗 (Bartlett window)


$$\begin {equation}\begin{aligned} \\w(n)= \left\{\begin{array}{ll} \frac{2n}{N-1}, & 0\leq n \leq (N-1)/2 \\ 2-\frac{2n}{N-1}, & N/2 \leq n \leq N-1 \\ 0,& n > N-1 \end{array}\right. \end{aligned} \label{eq:Bartlett_win_def} \end{equation}$$

布莱克曼窗 (Blackman window)


$$\begin {equation}\begin{aligned} \\w(n)= \left\{\begin{array}{ll} 0.42-0.5\cos(2\pi \frac{n}{N-1})+0.08\cos(4\pi \frac{n}{N-1}), & 0\leq n \leq N-1 \\ 0,& n > N-1 \end{array}\right. \end{aligned} \label{eq:Blackman_win_def} \end{equation}$$

实验代码

以下代码，绘制了 Hanning ，Hamming， Bartlett， Blackman 窗的时域波形和幅度频谱。

close all; clear; clc;
dtft_calc=@(xn_col_vec,w_delta) ... 
    exp(-j*(([0:w_delta:2*pi].')*[0:1:length(xn_col_vec)-1])) * xn_col_vec;

N = 40;
n_idx = [0:N-1].';
xn1 = hanning(N);  s1 = 'hanning  '; 
xn2 = hamming(N);  s2 = 'hamming  '; 
xn3 = blackman(N); s3 = 'blackman ';  
xn4 = bartlett(N); s4 = 'bartlett ';  

w_delta = 0.0001*pi;

yn1_dB = 20*log10(abs(fftshift(dtft_calc(xn1, w_delta)))); yn1_dB_norm  = yn1_dB - max(yn1_dB);
yn2_dB = 20*log10(abs(fftshift(dtft_calc(xn2, w_delta)))); yn2_dB_norm  = yn2_dB - max(yn2_dB);
yn3_dB = 20*log10(abs(fftshift(dtft_calc(xn3, w_delta)))); yn3_dB_norm  = yn3_dB - max(yn3_dB);
yn4_dB = 20*log10(abs(fftshift(dtft_calc(xn4, w_delta)))); yn4_dB_norm  = yn4_dB - max(yn4_dB);

w_idx = [-pi:w_delta:pi]/pi;
figure;  
subplot(4,2,[1 2 3 4]);plot(n_idx,xn1,':',n_idx,xn2,'--',n_idx,xn3,'-.', n_idx,xn4, ...
'lineWidth',1.5); grid on; legend({s1,s2,s3,s4}, 'FontSize',12);
title('Time Domain','FontSize',12);
subplot(4,2,5);plot(w_idx,yn1_dB_norm,'lineWidth',1.5);title([s1 'spectrum']); grid on;ylim([-160,0]); ylabel('Amplitude(dB)');xlabel('\omega: \pi rad');
subplot(4,2,6);plot(w_idx,yn2_dB_norm,'lineWidth',1.5);title([s2 'spectrum']); grid on;ylim([-160,0]); ylabel('Amplitude(dB)');xlabel('\omega: \pi rad');
subplot(4,2,7);plot(w_idx,yn3_dB_norm,'lineWidth',1.5);title([s3 'spectrum']); grid on;ylim([-160,0]); ylabel('Amplitude(dB)');xlabel('\omega: \pi rad');
subplot(4,2,8);plot(w_idx,yn4_dB_norm,'lineWidth',1.5);title([s4 'spectrum']); grid on;ylim([-160,0]); ylabel('Amplitude(dB)');xlabel('\omega: \pi rad');

实验作业—经典窗函数的频域分析

请使用以上代码，运行仿真，观察频谱波形，回答以下问题

凯泽窗函数

从前文的内容可以看到，各种经典窗函数之间的区别在于频响曲线的主瓣宽度和主瓣对旁瓣的抑制程度。 并且，从实验中可以看到，主瓣宽度和旁瓣抑制比是一对矛盾的参数，即：

另一方面，从信号分析的角度，此项调节参数的意义在于：

与之前的经典窗函数相比，凯泽窗（kaiser window）是一种更为先进的窗函数。凯泽窗拥有连续的形状调节能力， 它比之前的窗函数多了一个参数因子 β ， 通过改变 参数 β的数值，凯泽窗可以改变形状。 当 β = 0时，凯泽窗等价于 矩形窗函数，此时的主瓣宽度最窄，当增大 β值后，凯泽窗的 主瓣会展宽，但是可以获得更好的旁瓣抑制比。

凯泽窗的定义如下：


$$\begin {equation}\large\begin{aligned} \\w(n)= \left\{\begin{array}{ll} \frac{I_0\left[\beta\sqrt{1-(\frac{2n}{N-1}-1)^2}\right]}{I_0[\beta]}, & 0\leq n \leq N-1 \\ 0,& n > N-1 \end{array}\right. \end{aligned} \label{eq:Kaiser_win_def} \end{equation}$$

其中 I0[ ⋅ ] 表示 第一类零阶修正贝塞尔（Bessel）函数。这个函数的形式比较复杂，此处不展开讨论。

实验代码

以下代码用于计算给定N值下，多个β参数对应的凯泽窗函数的时域波形和分贝尺度的幅度谱。

close all; clear; clc;
dtft_calc=@(xn_col_vec,w_delta) ... 
    exp(-j*(([0:w_delta:2*pi].')*[0:1:length(xn_col_vec)-1])) * xn_col_vec;

N = 20;
n_idx = [0:N-1].';
xn1 = kaiser(N,0);  s1 = '\beta = 0 '; 
xn2 = kaiser(N,1);  s2 = '\beta = 2 '; 
xn3 = kaiser(N,2);  s3 = '\beta = 4 ';  
xn4 = kaiser(N,4);  s4 = '\beta = 16';  

w_delta = 0.0001*pi;

yn1_dB = 20*log10(abs(fftshift(dtft_calc(xn1, w_delta)))); yn1_dB_norm  = yn1_dB - max(yn1_dB);
yn2_dB = 20*log10(abs(fftshift(dtft_calc(xn2, w_delta)))); yn2_dB_norm  = yn2_dB - max(yn2_dB);
yn3_dB = 20*log10(abs(fftshift(dtft_calc(xn3, w_delta)))); yn3_dB_norm  = yn3_dB - max(yn3_dB);
yn4_dB = 20*log10(abs(fftshift(dtft_calc(xn4, w_delta)))); yn4_dB_norm  = yn4_dB - max(yn4_dB);

w_idx = [-pi:w_delta:pi]/pi;
figure;  
subplot(4,2,[1 2 3 4]);plot(n_idx,xn1,'-',n_idx,xn2,'--',n_idx,xn3,'-.', n_idx,xn4, ':', ...
'lineWidth',1.5); ylim([0,1.2]);grid on; legend({s1,s2,s3,s4}, 'FontSize',12);
title('Time Domain','FontSize',12);
subplot(4,2,5);plot(w_idx,yn1_dB_norm,'lineWidth',1.5);title([s1 'spectrum']); grid on;ylim([-100,0]); ylabel('Amplitude(dB)');xlabel('\omega: \pi rad');
subplot(4,2,6);plot(w_idx,yn2_dB_norm,'lineWidth',1.5);title([s2 'spectrum']); grid on;ylim([-100,0]); ylabel('Amplitude(dB)');xlabel('\omega: \pi rad');
subplot(4,2,7);plot(w_idx,yn3_dB_norm,'lineWidth',1.5);title([s3 'spectrum']); grid on;ylim([-100,0]); ylabel('Amplitude(dB)');xlabel('\omega: \pi rad');
subplot(4,2,8);plot(w_idx,yn4_dB_norm,'lineWidth',1.5);title([s4 'spectrum']); grid on;ylim([-100,0]); ylabel('Amplitude(dB)');xlabel('\omega: \pi rad');

实验作业——凯泽窗和经典窗的对比

短时傅里叶变换

　　时域和频域，是我们分析信号的两种手段。虽然离散傅里叶变换定义式的计算过程是在整个时间轴上进行， 但是现实世界中的信号往往是暂态信号，即，信号在一段时间内，其能量集中在一部分频率上。

　　短时傅里叶变换，给我们这样一种方法，它在一个无穷长的信号上，开启了一个有限长的时间窗口， 对于落入这个时间窗口的信号样本，将其作为一个有限长信号序列进行离散傅里叶变换。当时间窗口沿着时间轴逐个 样点的向后滑动过去，每滑动一个样点，则计算出一幅频谱曲线出来，于是，一维的时间信号就变成了二维 的时间/频率信号。

信号 x[n]的短时傅里叶变换的定义公式如下所示：


$$\begin {equation}\large\begin{aligned} X[n,\omega] = \sum_{m = -\infty}^{\infty}x[n+m]u[m]e^{-j\omega m} \end{aligned} \label{eq:st_dtft_def} \end{equation}$$

其中:

短时傅里叶变换，为我们提供了一个可以在时域和频域进行调节的信号分析方法：

实验代码 — 时频联合分析


$$\begin {equation}\large\begin{aligned} x[n] &= \sum_{i=0}^2 x_i[n] \\ x_i[n] &= \left\{\begin{array}{ll} A_i\sin(\omega_i(n-S_i)), & S_i\leq n \leq S_i+M_i-1 \\ 0,& n \geq S_i+M_i \end{array}\right. \end{aligned} \label{eq:x_x012} \end{equation}$$

close all; clear; clc;

dtft_calc=@(xn_col_vec,w_delta) ... 
    exp(-j*(([0:w_delta:2*pi].')*[0:1:length(xn_col_vec)-1])) * xn_col_vec;

win_len = 40; kaiser_beta = 10; 

N = 650; w_delta = 0.001*pi; 
w0 = 0.1*pi;  A0 =   1 ; M0 = 200; S0 = 100;
w1 = 0.2*pi;  A1 = 0.4 ; M1 = 150; S1 = 320; 
w2 = 0.5*pi;  A2 = 0.1 ; M2 = 100; S2 = 150; 

n_idx = [0:N-1].';

x0 = A0*sin([0:(M0 -1)].' * w0);
x1 = A1*sin([0:(M1 -1)].' * w1);
x2 = A2*sin([0:(M2 -1)].' * w2);

x0zN =zeros(N,1); x0zN(S0:(S0 + M0-1)) = x0;
x1zN =zeros(N,1); x1zN(S1:(S1 + M1-1)) = x1;
x2zN =zeros(N,1); x2zN(S2:(S2 + M2-1)) = x2;

xn = x0zN+ x1zN+ x2zN; xn = xn ./ max(abs(xn));
xn_win = xn .* kaiser(N,kaiser_beta);

yn = dtft_calc(xn_win, w_delta);
yn = fftshift(yn);
yn_dB = 20*log10(abs(yn));
yn_dB_norm  = yn_dB - max(yn_dB);
w_idx = [-pi:w_delta:pi]/pi;
figure;  
subplot(4,1,1 );plot(n_idx,xn  ); ylim([-1.1, 1.1]);grid on;
subplot(4,1,2 );plot(n_idx,x0zN); ylim([-1.1, 1.1]);grid on;
subplot(4,1,3 );plot(n_idx,x1zN); ylim([-1.0, 1.0]);grid on;
subplot(4,1,4 );plot(n_idx,x2zN); ylim([-1.0, 1.0]);grid on;
figure;
subplot(1,2,1);plot(w_idx,abs(yn)); grid on;
subplot(1,2,2);plot(w_idx,yn_dB_norm); grid on;ylim([-100,10]);

% calculate spectrogram 
ynspg = zeros(length(yn), N-win_len); 
xn_ii = zeros(win_len,1);
for(ii = 1:N) 
  xn_ii     = [xn_ii(2:win_len);xn(ii)];
  xn_ii_win = xn_ii .* kaiser(win_len,kaiser_beta);
  yn_ii = fftshift(dtft_calc(xn_ii_win, w_delta));
  yn_ii_dB = 20*log10(abs(yn_ii)+1E-10);
  ynspg(:,ii) = yn_ii_dB;
end
max_spg = max(max(ynspg));
ynspg_norm = ynspg - max_spg; % normalization spectrum
figure;
mesh([0:N-1],w_idx,ynspg_norm);
axis tight; view(0,90); colorbar;
xlabel('Time: n');ylabel('Frequency \omega: \pi rad');
title(['Spectrogram: win\_len is ',num2str(win_len), ', beta is ', num2str(kaiser_beta)]);

实验作业

使用上面的代码，运行仿真，设定仿真参数如下，调节 滑动窗长度 win_len 和 凯泽窗形状因子 beta，运行仿真，回答问题 
N = 650; w_delta = 0.001*pi; 
w0 = 0.1*pi;  A0 =   1 ; M0 = 200; S0 = 100;
w1 = 0.2*pi;  A1 = 0.4 ; M1 = 150; S1 = 320; 
w2 = 0.5*pi;  A2 = 0.1 ; M2 = 100; S2 = 150;

数字滤波器的幅频和相频响应

延迟叠加系统的频率响应—单频信号解析法研究

　　回顾《电路原理》 和 《信号与系统》课程，带有频率选择性的电路系统， 比如无源滤波器， 其组成部分中必须有储能元件，即：电感和电容， 只包含电阻 元件的电路网络是无法构成滤波器的。对于储能元件而言，其区别于电阻元件的最大 特征是具有电抗特性，即其输出信号相对于输入信号存在时间延迟。

延迟叠加系统对于单频信号的输出表达式

　　 对于数字信号而言，情况也类似，如果想构成带有频率选择性的滤波器，数字系统 中也必须带有对输入信号的延迟单元。 下面首先分析一个最简单的余弦信号的抽样 序列 x[n] = cos(ω0n) 通过一个延迟叠加系统的频率响应，即 ： y[n] = x[n] − x[n − D] 的频率响应。 首先，我们从时域来分析，


$$\begin {equation}\large\begin{aligned} x[n] &= \cos(\omega_0 n) \\ y[n] &= x[n] - x[n-D] \\ &= \cos(\omega_0 n)-\cos(\omega_0 (n-D)) \end{aligned} \end{equation}$$

根据三角函数的积化和差公式，有：
$$\begin {equation}\large\begin{aligned} \cos(\alpha)-\cos(\beta) &= -2\sin(\frac{\alpha+\beta}{2}) \sin(\frac{\alpha-\beta}{2}) \end{aligned} \end{equation}$$
则：
$$\begin {equation}\large\begin{aligned} y[n] &= \cos(\omega_0 n)-\cos(\omega_0 (n-D)) \\ &= -2\sin(\frac{2\omega_0 n-\omega_0 D}{2}) \sin(\frac{\omega_0 D}{2}) \\ &= \cos(\omega_0 n + \frac{\pi}{2}-\frac{\omega_0 D}{2} )\cdot 2 \sin(\frac{\omega_0 D}{2}) \end{aligned} \label{eq:cos_w_n_s_D} \end{equation}$$

延迟叠加系统的幅频和相频响应

观察公式 $\eqref{eq:cos_w_n_s_D}$ ，会发现:

频率—幅频、相频响应参数计算表

ω0 $ y[n]$ |H(ω)| Φ(ω)
0 0 0 0
Δ cos(n ⋅ Δ + π/2 − 2Δ) ⋅ 2sin(2Δ) |2sin(2Δ)| ≈ 4Δ π/2 − 2Δ ≈ π/2
π/8 cos(n ⋅ π/8 + π/2 − 2π/8) ⋅ 2sin(2π/8) $|2\sin(2\pi/8) |=\sqrt{2}=1.414$ π/2 − 2π/8 = 2π/8
2π/8 cos(n ⋅ 2π/8 + π/2 − 4π/8) ⋅ 2sin(4π/8) |2sin(4π/8)| = 2 π/2 − 4π/8 = 0
3π/8 cos(n ⋅ 3π/8 + π/2 − 6π/8) ⋅ 2sin(6π/8) $|2\sin(6\pi/8) |=\sqrt{2}=1.414$ π/2 − 6π/8 =  − 2π/8
4π/8 − Δ cos(n ⋅ 4π/8 + π/2 − 8π/8 + 2Δ) ⋅ 2sin(8π/8 − 2Δ) |2sin(8π/8 − 2Δ)| ≈ 4Δ π/2 − 8π/8 + 2Δ =  − 4π/8 + 2Δ ≈  − π/2
4π/8 cos(n ⋅ 4π/8 + π/2 − 8π/8) ⋅ 2sin(8π/8) |2sin(8π/8)| = 0 0
4π/8 + Δ cos(n ⋅ 4π/8 + π/2 − 8π/8 − 2Δ) ⋅ 2sin(8π/8 + 2Δ) |2sin(8π/8 + 2Δ)| ≈ 4Δ π/2 − 8π/8 − 2Δ + π = π/2 − 2Δ ≈ π/2
5π/8 cos(n ⋅ 5π/8 + π/2 − 10π/8) ⋅ 2sin(10π/8) $|2\sin(10\pi/8) |=\sqrt{2}=1.414$ π/2 − 10π/8 + π = 2π/8
6π/8 cos(n ⋅ 6π/8 + π/2 − 12π/8) ⋅ 2sin(12π/8) |2sin(12π/8)| = 2 π/2 − 12π/8 + π = 0
7π/8 cos(n ⋅ 7π/8 + π/2 − 14π/8) ⋅ 2sin(14π/8) $|2\sin(14\pi/8) |=\sqrt{2}=1.414$ π/2 − 14π/8 + π =  − 2π/8
8π/8 − Δ cos(n ⋅ 8π/8 + π/2 − 16π/8 + 2Δ) ⋅ 2sin(16π/8 − 2Δ) |2sin(16π/8 − 2Δ)| ≈ 4Δ π/2 − 16π/8 + 2Δ + π =  − 4π/8 + 2Δ

延迟叠加系统的零极点分析

延迟叠加系统的Z变换系统函数推导如下：
$$\begin {equation}\begin{aligned} y[n] &= x[n] - x[n - D] \\ Y(Z) &= X(Z)\cdot (1 - Z^{-D}) \\ H(Z) &= \frac{Z^D-1}{Z^D} \end{aligned} \label{eq:yn_xn_HZ_D} \end{equation}$$

根据复变函数的知识， 系统函数 H(Z) 有 D个 等分圆周的一阶零点，和 1 个 位于原点的D阶极点

延迟叠加系统的仿真验证

下面使用Matlab仿真代码来验算一下我们用的解析法得到的结论，以下代码用于绘制

为了增强曲线图的可读性，代码中使用了一些绘图命令

close all; clear; clc;
N   = 100; D   = 4; w0 = 1*pi/8; phy =0*pi/8; 
idx = [0:N-1].'; ZD = zeros(D-1,1);
xn  = cos(w0*idx + phy);
hn   = [1;ZD;-1];
xnD = cos(w0*(idx-D) + phy);
yn = xn - xnD;

figure;
plot(idx,xn ,'--', idx,xnD,':', idx,yn,'-','LineWidth',2); grid on;title('time signal');
legend('xn','xnD','yn');

[hw,w] = freqz(hn);

figure;
xs = strcat(num2str(([0:2*D].')), repmat(['\pi/' num2str(2*D)],2*D+1,1)); xss = size(xs);
xc = mat2cell(xs,ones(2*D+1,1),[xss(2)]); xc{1} = '0';
subplot(2,1,1); plot(w, abs(hw));grid on; ylim([0 2.1]);title('H(w) abs');
set(gca,'XTick',0:pi/(D*2):pi);set(gca,'XTickLabel', xc);
ys = strcat(num2str(([-D:D].')), repmat(['\pi/' num2str(2*D)],2*D+1,1)); yss = size(ys);
yc = mat2cell(ys,ones(2*D+1,1),[yss(2)]);yc{D+1} = '0';
subplot(2,1,2); plot(w, angle(hw)); grid on;title('H(w) phase')
set(gca,'XTick',0:pi/(D*2):pi);set(gca,'YTick',-pi/2:pi/(D*2):pi/2);
set(gca,'XTickLabel', xc);set(gca,'YTickLabel', yc); 
figure;
subplot(2,1,1); plot(w, real(hw));grid on; title('H(w) real');
set(gca,'XTick',0:pi/(D*2):pi); set(gca,'XTickLabel', xc);
subplot(2,1,2); plot(w, imag(hw));grid on; title('H(w) imag');
set(gca,'XTick',0:pi/(D*2):pi);
set(gca,'TickLabelInterpreter', 'tex'); set(gca,'XTickLabel', xc);
figure; zplane(roots(hn), roots([1;zeros(D,1)]));title('Z-plane');

仿真结果的幅频-相频分析

以上的仿真 给出了系统的幅频-相频曲线，以及给出了延迟叠加系统的输入、输出的信号波形。 接下来我们分析该系统的幅频和相频特征。

实验作业

线性系统的群延迟分析

群延迟的定义和物理意义

群延迟(group delay) 是用来表征系统的相频曲线的线性度的参数。 设系统函数 H(ω) 的 相频函数为 Φ(ω) ， 则 H(ω)群延时函数 grd [H(ω)] 定义如下：


$$\begin {equation}\begin{aligned} grd~[H(\omega)] = -\frac{d\Phi(\omega)}{d\omega} \end{aligned} \label{eq:grd_def} \end{equation}$$

即：群延迟函数是系统的相频函数对频率变量求导的负数，加上负号的原因是， 因果系统的输出总是滞后于输入， 则相频曲线总是呈现出负斜率， 加上负号后，群延迟则变成一个正值函数，便于我们分析观察。

正如其名字所示， 群延迟（group delay），是一种从时间延迟的角度描述系统的相移特征的方法。 设系统在 ω0 的相移为 Φ(ω0)，以及在 ω1 的相移为 Φ(ω1) 以及系统的群延迟频率函数为 grd(ω)，则有：


$$\begin {equation}\begin{aligned} \Phi(\omega_1) = \Phi(\omega_0) - \int^{\omega_1}_{\omega_0} grd(\omega) d\omega \end{aligned} \label{eq:grd_w0_w1} \end{equation}$$

对于群延迟函数，有一种情况是我们重点关注的， 即线性相位系统，Φ(ω) 是一条斜率为负数的直线， 此时 grd(ω) 函数是一个常数 τ,  (τ > 0) 则式 $\eqref{eq:grd_w0_w1}$变为：


$$\begin {equation}\begin{aligned} \Phi(\omega_1) = \Phi(\omega_0) - (\omega_1 - \omega_0)\cdot \tau \end{aligned} \label{eq:grd_w0_w1_tau} \end{equation}$$

则，当频率为 ω1 的余弦信号 x[n] = cos(ω1n) 通过群延迟为 τ,  (τ > 0) 的线性系统 H(ω) 时， 其输出信号 y[n] 可以表示为：


$$\begin {equation}\begin{aligned} \\ y[n] &= |H(\omega_1)|\cdot \cos(\omega_1 n + \Phi(\omega_1)) \\ &= |H(\omega_1)|\cdot \cos(\omega_1 n + \Phi(\omega_0) - (\omega_1 - \omega_0)\cdot \tau) \\ &= |H(\omega_1)|\cdot \cos(\omega_1 n - \omega_1\tau + \Phi(\omega_0) + \omega_0\tau) \\ &= |H(\omega_1)|\cdot \cos(\omega_1(n - \tau) + \Phi(\omega_0) + \omega_0\tau) \end{aligned} \label{eq:grd_y_omega_1_tau} \end{equation}$$

令 β = Φ(ω0) + ω0τ ， 则式 $\eqref{eq:grd_y_omega_1_tau}$ 变为:
$$\begin {equation}\begin{aligned} \\y[n] = |H(\omega_1)|\cdot \cos(\omega_1(n - \tau) + \beta) \end{aligned} \label{eq:grd_y_omega_1_tau_beta} \end{equation}$$

观察 式. $\eqref{eq:grd_y_omega_1_tau_beta}$， 可见：

群延迟与频域分量的时间色散

线性相位系统有一个重要特性，即：不会对输入信号引发时间色散。详见以下分析：

广义线性相位系统

回顾之前延迟叠加系统仿真实验， 用 freqz 函数绘制的延迟叠加系统的相频函数曲线，以及我们之前推导的相频函数曲线。 会发现在每一个幅频响应为0的位置，相频函数会发生跳变，这就是造成相频函数曲线不连续的原因。但是在 两个跳变之间的相频函数仍然是一段直线，因此，延迟叠加系统是一个广义线性相位系统。

下面我们使用仿真实验来观察广义线性相位系统对个频率分量的时间色散特性。

仿真实验：线性相位与时间色散

代码说明

实验代码 
close all; clc; clear;

N   = 300; D   = 20; w0 = 2*pi/100; w1 = 3*pi/100; w2 = 15*pi/100;
w012 = [w0;w1;w2]; idx = [0:N-1].'; ZD = zeros(D-1,1);
hn   = [1;ZD;-1];
[hw,w] = freqz(hn,1,1000);
fn = w/pi;
ang_hw = angle(hw)/pi;

set(groot, 'DefaultLineLineWidth',2);
set(groot, 'DefaultAxesFontSize', 14);
set(groot, 'DefaultAxesXGrid', 'on');
set(groot, 'DefaultAxesYGrid', 'on');
set(groot, 'DefaultAxesXMinorGrid', 'on');
set(groot, 'DefaultAxesYMinorGrid', 'on');

figure;
subplot(2,1,1);plot(fn,abs(hw)) ;grid on; 
xlabel('\omega (\pi rad)'); ylabel('|H(\omega)|');  title('Amp-Freq');
subplot(2,1,2);plot(fn, ang_hw);grid on;
xlabel('\omega (\pi rad)'); ylabel('\Phi(\omega) (\pi rad)');title('Phase-Freq');

xn_c_012  = zeros(N,3); xnD_c_012 = zeros(N,3); yn_c_012  = zeros(N,3); 
xn_s_012  = zeros(N,3); xnD_s_012 = zeros(N,3); yn_s_012  = zeros(N,3);
ang_xn_012 = zeros(N,3);ang_xnD_012 = zeros(N,3);ang_yn_012 = zeros(N,3);

for(ii = 1:3)
  xn_c_012(:,ii) = cos(w012(ii)*idx); xnD_c_012(:,ii)= cos(w012(ii)*(idx-D) ); 
  xn_s_012(:,ii) = sin(w012(ii)*idx); xnD_s_012(:,ii)= sin(w012(ii)*(idx-D) );
  yn_c_012(:,ii) = xn_c_012(:,ii) - xnD_c_012(:,ii); yn_c_012(:,ii) = yn_c_012(:,ii)/max(abs(yn_c_012(:,ii)));
  yn_s_012(:,ii) = xn_s_012(:,ii) - xnD_s_012(:,ii); yn_s_012(:,ii) = yn_s_012(:,ii)/max(abs(yn_s_012(:,ii)));
  ang_xn_012 (:,ii) = (angle(xn_c_012 (:,ii)+j*xn_s_012 (:,ii)))/pi; 
  ang_xnD_012(:,ii) = (angle(xnD_c_012(:,ii)+j*xnD_s_012(:,ii)))/pi; 
  ang_yn_012 (:,ii) = (angle(yn_c_012 (:,ii)+j*yn_s_012 (:,ii)))/pi; 
end

x1c=  xn_c_012(:,1); x1c_L = 'xn_c_0'; x1s=  xn_s_012(:,1);  x1s_L =  'xn_s_0';
x2c=  xn_c_012(:,2); x2c_L = 'xn_c_1'; x2s=  xn_s_012(:,2);  x2s_L =  'xn_s_1';
x3c=  xn_c_012(:,3); x3c_L = 'xn_c_2'; x3s=  xn_s_012(:,3);  x3s_L =  'xn_s_2';
d1c= xnD_c_012(:,1); d1c_L = 'xnD_c_0';d1s= xnD_s_012(:,1);  d1s_L = 'xnD_s_0';
d2c= xnD_c_012(:,2); d2c_L = 'xnD_c_1';d2s= xnD_s_012(:,2);  d2s_L = 'xnD_s_1';
d3c= xnD_c_012(:,3); d3c_L = 'xnD_c_2';d3s= xnD_s_012(:,3);  d3s_L = 'xnD_s_2';
y1c=  yn_c_012(:,1); y1c_L =  'yn_c_0';y1s=  yn_s_012(:,1);  y1s_L =  'yn_s_0';
y2c=  yn_c_012(:,2); y2c_L =  'yn_c_1';y2s=  yn_s_012(:,2);  y2s_L =  'yn_s_1';
y3c=  yn_c_012(:,3); y3c_L =  'yn_c_2';y3s=  yn_s_012(:,3);  y3s_L =  'yn_s_2';
x1a=  ang_xn_012(:,1); x1a_L =  'Ang_xn_0';d1a= ang_xnD_012(:,1); d1a_L = 'Ang_xnD_0';
x2a=  ang_xn_012(:,2); x2a_L =  'Ang_xn_1';d2a= ang_xnD_012(:,2); d2a_L = 'Ang_xnD_1';
x3a=  ang_xn_012(:,3); x3a_L =  'Ang_xn_2';d3a= ang_xnD_012(:,3); d3a_L = 'Ang_xnD_2';
y1a=  ang_yn_012(:,1); y1a_L =  'Ang_yn_0';
y2a=  ang_yn_012(:,2); y2a_L =  'Ang_yn_1';
y3a=  ang_yn_012(:,3); y3a_L =  'Ang_yn_2';

txc =  'xn cos signal';txs =  'xn sin signal'; txa = 'xn  Angle';
tdc = 'xnD cos signal';tds = 'xnD sin signal'; tda = 'xnD Angle';
tyc =  'yn cos signal';tys =  'yn sin signal'; tya = 'yn  Angle';

figure;
subplot(3,1,1); plot(idx,x1c, idx,x2c, idx,x3c); title(txc) ; legend(x1c_L,x2c_L,x3c_L);
subplot(3,1,2); plot(idx,d1c, idx,d2c, idx,d3c); title(tdc) ; legend(d1c_L,d2c_L,d3c_L);
subplot(3,1,3); plot(idx,y1c, idx,y2c, idx,y3c); title(tyc) ; legend(y1c_L,y2c_L,y3c_L);
figure;
subplot(3,1,1); plot(idx,x1s, idx,x2s, idx,x3s); title(txs) ; legend(x1s_L,x2s_L,x3s_L);
subplot(3,1,2); plot(idx,d1s, idx,d2s, idx,d3s); title(tds) ; legend(d1s_L,d2s_L,d3s_L);
subplot(3,1,3); plot(idx,y1s, idx,y2s, idx,y3s); title(tys) ; legend(y1s_L,y2s_L,y3s_L);
figure; 
yLs = '\Phi(n) (\pi rad)';
subplot(3,1,1); plot(idx,x1a, idx,x2a, idx,x3a); title(txa ) ; legend(x1a_L,x2a_L,x3a_L); xlabel('n'); ylabel(yLs);
subplot(3,1,2); plot(idx,d1a, idx,d2a, idx,d3a); title(tda ) ; legend(x1a_L,x2a_L,x3a_L); xlabel('n'); ylabel(yLs);
subplot(3,1,3); plot(idx,y1a, idx,y2a, idx,y3a); title(tya ) ; legend(x1a_L,x2a_L,x3a_L); xlabel('n'); ylabel(yLs);

实验作业

FIR 数字滤波器

MATLAB的三种常用滤波器设计函数

Matlab 提供了三种 FIR 滤波器设计方法

窗函数法设计FIR及实验

clear;clc;close all
set(groot, 'DefaultLineLineWidth',1.5); set(groot, 'DefaultAxesFontSize', 14); 

N_tap = 32; f_cutoff = 0.6; beta1 = 0;beta2 = 5;

N_freqz_pt = 1000;
b1 = fir1(N_tap-1,f_cutoff,kaiser(N_tap,beta1));
[h1,w1]  = freqz(b1,1,N_freqz_pt);
h1_abs = abs(h1);
h1_dB = 20*log10(h1_abs+1E-10);

b2 = fir1(N_tap-1,f_cutoff,kaiser(N_tap,beta2));
[h2,w2]  = freqz(b2,1,N_freqz_pt);
h2_abs = abs(h2);
h2_dB = 20*log10(h2_abs+1E-10);
ttstr1 = [num2str(N_tap), 'Taps, beta = ' , num2str(beta1), ', Linear Scale'];
ttstr2 = [num2str(N_tap), 'Taps, beta = ' , num2str(beta1), ', dB Scale'    ];
ttstr3 = [num2str(N_tap), 'Taps, beta = ' , num2str(beta2), ', Linear Scale'];
ttstr4 = [num2str(N_tap), 'Taps, beta = ' , num2str(beta2), ', dB Scale'    ];
figure;
subplot(2,2,1);plot(w1/pi,h1_abs); grid on; title(ttstr1);ylabel('Amplitude ');xlabel('\omega / \pi ');
subplot(2,2,3);plot(w1/pi,h1_dB ); grid on; title(ttstr2);ylabel('Amplitude ');xlabel('\omega / \pi ');
subplot(2,2,2);plot(w1/pi,h2_abs); grid on; title(ttstr3);ylabel('Amplitude ');xlabel('\omega / \pi ');
subplot(2,2,4);plot(w1/pi,h2_dB ); grid on; title(ttstr4);ylabel('Amplitude ');xlabel('\omega / \pi ');

figure;
idx = [0:N_tap-1];
subplot(2,1,1);stem(idx,b1,'o','filled','linestyle','none');grid on;hold; plot(idx,b1);
subplot(2,1,2);stem(idx,b2,'o','filled','linestyle','none');grid on;hold; plot(idx,b2);

fvtool(b1,1); fvtool(b2,1);

实验作业

请尝试不同的窗函数 beta 因子，用最少的抽头数量达到要求

频率抽样法设计FIR及实验

clear;clc;close all
set(groot, 'DefaultLineLineWidth',1.5);
set(groot, 'DefaultAxesFontSize', 14);

N_tap1 =  32; 
N_tap2 = 128; 
N_freqz_pt = 1000*N_tap1;

f1 = [0 0.1 ]; f2 = [0.1:0.001:0.5]   ; f3 = [0.501  0.6 0.7 0.8 0.9 1];
m1 =  [0 0  ]; m2 = f2*1/(0.5-0.1)-0.1; m3 = [0.0    0.0 0.0 0.0 0.0 0];

m2_dB = 20*log10(m2); 

f    = [f1 f2    f3];
m    = [m1 m2    m3];
m_dB = [m1 m2_dB m3];

b1 = fir2(N_tap1-1,f,m);
b2 = fir2(N_tap2-1,f,m);
[h1,w1]  = freqz(b1,1,N_freqz_pt);
[h2,w2]  = freqz(b2,1,N_freqz_pt);
h1_abs = abs(h1);
h2_abs = abs(h2);
h1_dB = 20*log10(h1_abs+1E-10);
h2_dB = 20*log10(h2_abs+1E-10);

ttstr1 = [ 'Filter frequency reponse ' num2str(N_tap1), ' Taps '];
ttstr2 = [ 'Filter frequency reponse ' num2str(N_tap2), ' Taps '];
ystr = 'Amplitude, Linear Scale';
figure;
subplot(2,1,1);plot(f,m,'--', w1/pi,h1_abs); grid on; title([ttstr1,', Linear']);ylabel(ystr);xlabel('\omega / \pi ');legend('Ideal', 'Derived');
subplot(2,1,2);plot(f,m,'--', w1/pi,h2_abs); grid on; title([ttstr2,', Linear']);ylabel(ystr);xlabel('\omega / \pi ');legend('Ideal', 'Derived');
figure;
subplot(2,1,1);plot(f,m_dB,'--', w1/pi,h1_dB); grid on; title([ttstr1,', dB']);ylabel(ystr);xlabel('\omega / \pi ');legend('Ideal', 'Derived');
subplot(2,1,2);plot(f,m_dB,'--', w1/pi,h2_dB); grid on; title([ttstr2,', dB']);ylabel(ystr);xlabel('\omega / \pi ');legend('Ideal', 'Derived');

figure;
subplot(2,1,1);plot([0:N_tap1-1], b1);grid;title( [num2str(N_tap1), ' Taps ']);
subplot(2,1,2);plot([0:N_tap2-1], b2);grid;title( [num2str(N_tap2), ' Taps ']);

实验作业

Parks-McClellan法 设计最优化FIR及实验

实验代码


clear;clc;close all
set(groot, 'DefaultLineLineWidth',1.5);
set(groot, 'DefaultAxesFontSize', 14);


N_tap1 =  30; 
N_tap2 =  60; 
N_freqz_pt = 1000*N_tap1;


f = [0   0.1  0.15  1];
m = [1.0 1.0  0.0   0];

b1 = firpm(N_tap1-1,f,m);
b2 = firpm(N_tap2-1,f,m);
[h1,w1]  = freqz(b1,1,N_freqz_pt);
[h2,w2]  = freqz(b2,1,N_freqz_pt);
h1_abs = abs(h1);
h2_abs = abs(h2);
h1_dB = 20*log10(h1_abs+1E-10);
h2_dB = 20*log10(h2_abs+1E-10);

ttstr1 = [ 'Filter frequency reponse ' num2str(N_tap1), ' Taps '];
ttstr2 = [ 'Filter frequency reponse ' num2str(N_tap2), ' Taps '];
ystr = 'Amplitude, Linear Scale';
figure;
subplot(2,1,1);plot(f,m,'--', w1/pi,h1_abs); grid ; title(ttstr1);ylabel(ystr);xlabel('\omega / \pi ');legend('Ideal', 'Derived');
subplot(2,1,2);plot(f,m,'--', w1/pi,h2_abs); grid ; title(ttstr2);ylabel(ystr);xlabel('\omega / \pi ');legend('Ideal', 'Derived');

figure;
subplot(2,1,1);plot([0:N_tap1-1], b1);grid; hold on; stem([0:N_tap1-1], b1);
subplot(2,1,2);plot([0:N_tap2-1], b2);grid; hold on; stem([0:N_tap2-1], b2);

figure;
subplot(2,1,1);plot(w1/pi, h1_dB);grid 
subplot(2,1,2);plot(w1/pi, h2_dB);grid 
fvtool(b1,1);fvtool(b2,1);

实验作业

使用FIR滤波器进行信号滤波

Matlab的filter函数

滤波实验代码及作业

clear;clc;close all
set(groot, 'DefaultLineLineWidth',1.5); set(groot, 'DefaultAxesFontSize', 14); 


N_tap = 64; f_cutoff = 0.1; beta1 = 0;;

N_freqz_pt = 1000;

b1 = fir1(N_tap-1,f_cutoff,kaiser(N_tap,beta1));
[h1,w1]  = freqz(b1,1,N_freqz_pt);
h1_abs = abs(h1);
h1_dB = 20*log10(h1_abs+1E-10);

ttstr1 = [num2str(N_tap), 'Taps, beta = ' , num2str(beta1), ', Linear Scale'];
ttstr2 = [num2str(N_tap), 'Taps, beta = ' , num2str(beta1), ', dB Scale'    ];
figure;
idx = [0:N_tap-1];
subplot(3,1,1);plot(w1/pi,h1_abs); grid on; title(ttstr1);ylabel('Amplitude ');xlabel('\omega / \pi ');
subplot(3,1,2);plot(w1/pi,h1_dB ); grid on; title(ttstr2);ylabel('Amplitude ');xlabel('\omega / \pi ');
subplot(3,1,3);stem(idx,b1,'o','filled','linestyle','none');grid on;hold; plot(idx,b1);title('b1');

fvtool(b1,1); 

Nx = 300; w1 = 0.05*pi; w2 = 0.1*pi ;
idx_Nx = [0:Nx-1].';
x1 = sin(w1* idx_Nx);
x2 = sin(w2* idx_Nx);

y1 = filter(b1,1,x1);
y2 = filter(b1,1,x2);

figure;
subplot(2,1,1);plot(idx_Nx,x1,idx_Nx,y1);grid;legend('x1','y1');
subplot(2,1,2);plot(idx_Nx,x2,idx_Nx,y2);grid;legend('x2','y2');

实验作业