数字信号处理-基础一

数字频率与模拟频率的关系与特点

link1link2link3

模拟频率、数字角频率和DFT频率分辨率

数字频率与模拟频率的定义

模拟频率ff:每秒经历多少个周期,单位为HzHz

模拟角频率Ω\Omega:每秒经历多少弧度,单位为rad/srad/s

数字角频率ω\omega:采样点之间的弧度,单位为radrad

数字信号是由模拟信号采样而来,采样频率不一样,采样点的时间就不一样。因此用每秒经历多少个周期已无多大意义,所以。

{fs1Ω=2πfrad/sω=Ω/fs=ΩTsradω=Ω/fs=2πf/fs(11)\begin{cases} f{\to}s^{-1}\\ \Omega=2{\pi}f{\to}rad/s\\ \omega={\Omega}/f_s={\Omega}T_s{\to}rad\\ \end{cases} \\ \omega={\Omega}/f_s=2{\pi}{\color{blue}f/fs}{\qquad}(1-1)

对于上式(1-1)的解释:

  • 数字角频率ω\omega是模拟角频率Ω\Omega对采样频率fsf_s的归一化;
  • f/fsf/f_s是一个无量纲的数,2π2{\pi}代表着弧度;
  • 即此频率(数字角频率:radrad)非彼频率(模拟角频率:rad/srad/s)。

数字角频率与采样频率有关:

(...,2π,π,0,π,2π,...)(...,fs,0.5fs,0,0.5fs,fs,...)(...,-2\pi,-\pi,0,\pi,2\pi,...)\\(...,fs,-0.5fs,0,0.5fs,fs,...)

为什么模拟角频率和数字角频率不一样

一个单位是rad/srad/s,另一个是radrad,肯定就不一样了啊。

模拟角频率Ω(,+)\Omega{\subset}(-\infin,+\infin),而数字角频率ω(π,π)\omega{\subset}(-\pi,\pi),当然也可以是(0,2π)(0,2\pi)。但由于数字角频率是具有周期性的,所以也可以认为数字角频率ω(,+)\omega{\subset(-\infin,+\infin)},只不过是周期性的。

fs=1Hzf_s=1Hz,当Ω=π/8\Omega=\pi/8Ω=17π/8\Omega=17\pi/8时,抽样序列如下:可以看到虽然模拟角频率增加了2π2\pi,但是由于采样点数和采样值都相同,所以实际的离散序列也是一样的。这也体现出了离散序列的角频率是以2π2\pi为周期的。

MATLAB代码如下:

step = 64;  % 用于产生模拟信号的精度
t = 0:1/(2*step):20;
t = t';
w1 = pi / 8;
x1 = cos(w1 .* t);
w2 = pi * 17 / 8;
x2 = cos(w2 .* t);

fs = 1;
ts = 0:1/fs:20;
ts = ts';
% 下面的作用是查找ts元素离t中最近元素的索引
D = abs(bsxfun(@minus, ts.', t));
M = min(D, [], 1);
[Index, ~] = find(bsxfun(@eq, M, D));
x1n = x1(Index);  % x1的抽样序列
x2n = x2(Index);  % x2的抽样序列

subplot(211); plot(t, x1, 'b'); 
hold on; stem(ts, x1n, 'r');
title('$$\Omega=\frac{\pi}{8}$$', 'Interpreter', 'latex');
subplot(212); plot(t, x2, 'b'); 
hold on; stem(ts, x2n, 'r');
title('$$\Omega=\frac{17\pi}{8}$$', 'Interpreter', 'latex');
xlabel(['fs=', num2str(fs), 'Hz']);

时域离散信号和时域离散系统

信号、系统

时域离散信号

常用的典型序列

MATLAB产生各种典型序列

序列的运算

加、减、乘、除、反转等。

时域离散系统

线性系统

当系统T的输入为x1(n)x_1(n)时,输出是y1(n)y_1(n);输入为x2(n)x_2(n)时,输出是y2(n)y_2(n)。若满足线性,则输入为x1(n)+x2(n)x_1(n)+x_2(n)时,输出应为y1(n)+y2(n)y_1(n)+y_2(n);且输入为ax1(n)ax_1(n)时,输出是ay1(n)ay_1(n)

if{y1(n)=T[x1(n)]y2(n)=T[x2(n)]have{T[x1(n)+x2(n)]=y1(n)+y2(n)T[ax1(n)]=ay1(n)if{\quad} \begin{cases} y_1(n)=T\big[x_1(n)\big]\\ y_2(n)=T\big[x_2(n)\big]\\ \end{cases} \\ have{\quad} \begin{cases} T\big[x_1(n)+x_2(n)\big]=y_1(n)+y_2(n)\\ T\big[ax_1(n)\big]=ay_1(n) \end{cases}

时不变系统

ify(n)=T[x(n)]havey(nn0)=T[x(nn0)]if{\quad}y(n)=T\big[x(n)\big]\\ have{\quad}y(n-n_0)=T\big[x(n-n_0)\big]

线性时不变系统及其输入与输出之间的关系

h(n)=T[δ(n)]x(n)=m=+x(m)δ(nm)y(n)=T[m=+x(m)δ(nm)]=m=+x(m)[δ(nm)]y(n)=m=+x(m)h(nm)=x(n)h(n)系统单位冲激响应:{\quad}h(n)=T\Big[\delta(n)\Big]\\ 输入信号:x(n)=\sum_{m=-\infin}^{+\infin}x(m)\delta(n-m)\\ 则输出为:y(n)=T\Big[\sum_{m=-\infin}^{+\infin}x(m)\delta(n-m)\Big]\\ =\sum_{m=-\infin}^{+\infin}x(m)\Big[\delta(n-m)\Big]\\ 根据时不变性质有:\\ y(n)=\sum_{m=-\infin}^{+\infin}x(m)h(n-m)=x(n)*h(n)

系统的因果性和稳定性

系统当前的输出只与之前的输入有关,后之后的输入无关。

h(n)=0n<0h(n)=0{\quad}n<0

时域离散系统的输入输出描述法—线性常系数差分方程

线性常系数差分方程

y(n)=i=0Mbix(ni)i=1Naiy(ni)ori=0Naiy(ni)=i=0Mbix(ni)a0=1y(n)=\sum_{i=0}^{M}b_ix(n-i)-\sum_{i=1}^{N}a_iy(n-i)\\ or{\quad}\sum_{i=0}^{N}a_iy(n-i)=\sum_{i=0}^{M}b_ix(n-i){\quad}a_0=1

差分方程的阶数是用方程y(ni)y(n-i)项中ii最大取值与最小取值之差确定的,即上式是N阶线性常系数差分方程。

线性常系数差分方程的求解

  • 经典解法:齐次解、特解
  • 递推解法
  • 交换域方法:将差分方程变换到zz
  • 还可以由差分方程求出系统的单位脉冲响应,然后与已知的输入序列进行卷积运算。

模拟信号数字处理方法

采样定理及A/D变换

xa(t)x_a(t)是模拟信号、pδ(t)p_{\delta}(t)是单位冲激序列、xa^(t)\hat{x_a}(t)是经理想采样后的信号。

pδ(t)=n=+δ(tnT)xa^(t)=xa(t)pδ(t){Xa(jω)=F[xa(t)]Xa^(jω)=F[xa^(t)]Pδ(jω)=F[pδ(t)]p_{\delta}(t)=\sum_{n=-\infin}^{+\infin}\delta(t-nT)\\ \hat{x_a}(t)=x_a(t)p_{\delta}(t)\\ \begin{cases} X_a(j\omega)=\mathscr{F}\Big[x_a(t)\Big]\\ \hat{X_a}(j\omega)=\mathscr{F}\Big[\hat{x_a}(t)\Big]\\ P_{\delta}(j\omega)=\mathscr{F}\Big[p_{\delta}(t)\Big]\\ \end{cases}

下面推导理想采样信号的频谱与原模拟信号频谱的关系:

Xa^(jΩ)=12πXa(jΩ)Pδ(jΩ)...Xa^(jΩ)=1Tk=+Xa(jΩjkΩs)\hat{X_a}(j\Omega)=\frac{1}{2\pi}X_a(j\Omega)*P_{\delta}(j\Omega)\\ ...\\ \hat{X_a}(j\Omega)=\frac{1}{T}\sum_{k=-\infin}^{+\infin}X_a(j\Omega-jk\Omega_s)

上式表明理想采样信号的频谱原模拟信号频谱沿频率轴,每间隔采样角频率Ωs=2πFS=2πT\Omega_s=2{\pi}F_S=\frac{2{\pi}}{T}重复出现一次,并叠加而形成的周期函数。或者说理想采样信号的频谱是原模拟信号的频谱以Ωs\Omega_s为周期,进行周期性延拓而成的。

因此,若要理想采样信号的频谱不重合,则需要满足:Ωswmax\Omega_s{\ge}w_{max}FsfmaxF_s{\ge}f_{max},这就是采样定理

数字序列转换成模拟信号

如果x(n)x(n)是在满足抽样定理的条件下得到的,那么可以通过一个低通滤波器不失真地将原模拟信号xa(t)x_a(t)恢复出来,低通滤波器的传输函数如下:

G(jΩ)={TΩ<12Ωs0Ω  12ΩsG(j\Omega)= \begin{cases} T & \big|\Omega\big|<\frac{1}{2}\Omega_s\\ 0 & \big|\Omega\big|\ {\ge}\ \frac{1}{2}\Omega_s\\ \end{cases}

其时域表达式为:

g(t)=12π+G(jΩ)ejΩtdΩ=12πΩs/2+Ωs/2TejΩtdΩ=sin(Ωst/2)Ωst/2Ωs=2πFs=2π/Tg(t)=sin(πt)/Tπt/Tg(t)=\frac{1}{2\pi}\int_{-\infin}^{+\infin}G(j\Omega)e^{j{\Omega}t}d_{\Omega}=\frac{1}{2\pi}\int_{-\Omega_s/2}^{+\Omega_s/2}Te^{j{\Omega}t}d_{\Omega}=\frac{\sin({\Omega_s}t/2)}{{\Omega_s}t/2}\\ {\because}{\qquad}{\Omega_s}=2{\pi}F_s=2{\pi}/T\\ {\therefore}g(t)=\frac{\sin({\pi}t)/T}{{\pi}t/T}

g(t)g(t)被称为内插函数。

x(n)=xa(nT)x(n)=x_a(nT)是经xa(t)x_a(t)满足采样定理采样后的离散信号,我们现在知道离散信号,如何恢复出模拟信号呢?转换公式如下(详细推导见DSP书第26页):

xa(t)=n=+xa(nT)sin[π(tnT)]/Tπ(tnT)/Tx_a(t)=\sum_{n=-\infin}^{+\infin}x_a(nT)\frac{\sin\Big[{\pi}(t-nT)\Big]\Big/T}{{\pi}(t-nT)\Big/T}

实际中采用D/AC完成数字信号到模拟信号的转换,包括三部分:解码器、零阶保持器和平滑滤波器。框图如下:

由时域离散信号xa(nT)x_a(nT)恢复模拟信号的过程是内插过程。

时域离散信号和系统的频域分析

时域离散信号的傅里叶变换的定义及性质

x(t)x(t)的傅里叶变换是X(ejω)X(e^{j\omega})x(n)x(n)的傅里叶变换也是X(ejω)X(e^{j\omega}),在频域都是连续的。

DTFT的定义

推导:

x(n)=x(t)δT(t)=x(t)n=δ(tnT)X(jω)=F[x(n)]=+x(t)n=δ(tnT)ejωtdt=n=+x(t)δ(tnT)ejωtdt=n=x(nT)ejωnT+δ(tnT)dt=n=x(nT)ejωnT{\because}{\quad}x(n)=x(t){\delta_T}(t)=x(t)\sum_{n=-\infin}^{\infin}{\delta}(t-nT)\\ {\therefore}{\quad} X(j\omega)={\mathscr{F}}[x(n)]=\int_{-\infin}^{+\infin}x(t)\sum_{n=-\infin}^{\infin}{\delta}(t-nT)e^{-j{\omega}t}d_t\\ =\sum_{n=-\infin}^{\infin}\int_{-\infin}^{+\infin}x(t)\delta(t-nT)e^{-j{\omega}t}d_t\\ =\sum_{n=-\infin}^{\infin}x(nT)e^{-j{\omega}nT}\int_{-\infin}^{+\infin}\delta(t-nT)d_t\\ =\sum_{n=-\infin}^{\infin}x(nT)e^{-j{\omega}nT}

T=1T=1,则变为:

x(n)=xa(nT)=x(t)t=nT{X(ejω)=n=+x(n)ejωnx(n)=12πππX(ejω)ejωndωx(n)=x_a(nT)=x(t)\Big|_{t=nT}\\ \begin{cases} X(e^{j\omega})=\sum_{n=-\infin}^{+\infin}x(n)e^{-j{\omega}n}\\ \\ x(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})e^{j{\omega}n}d_{\omega} \end{cases}

DTFT的性质

  • 周期性:X(ejω)=X(ejω+2πM)X(e^{j\omega})=X(e^{j{\omega+2{\pi}M}}),周期为2π2\pi

  • 线性

  • 时移和频移:

    [x(nn0)]=ejωn0X(ejω)[ejΩ0nx(n)]=X(ej(ωΩ0)) \Big[x(n-n_0)\Big]=e^{-j{\omega}n_0}X(e^{j\omega})\\ \Big[e^{j{\Omega_0}n}x(n)\Big]=X(e^{j({\omega}-{\Omega_0})})

  • 对称性:xe(n)x_e(n)为共轭对称序列(其实部是偶函数、虚部是奇函数,即满足xe(n)=xe(n)x_e(n)=x_e^{*}(-n));xo(n)x_o(n)为共轭反对称序列(其实部是奇函数、虚部是偶函数,即满足xo(n)=xo(n)x_o(n)=-x_o^{*}(-n))。

    x(n)=xr(n)+jxi(n)X(ejω)=Xe(ejω)+Xo(ejω)x(n)=xe(n)+xo(n)X(ejω)=XR(ejω)+jXI(ejω) x(n)={\color{blue}x_r(n)}+{\color{red}jx_i(n)}\\ X(e^{j\omega})={\color{blue}X_e(e^{j\omega})}+{\color{red}X_o(e^{j\omega})} \\ x(n)={\color{blue}x_e(n)}+{\color{red}x_o(n)}\\ X(e^{j\omega})={\color{blue}X_R(e^{j\omega})}+{\color{red}jX_I(e^{j\omega})}

  • 频域卷积定理:

y(n)=h(n)x(n)Y(ejω)=12πH(ejω)X(ejω)=12πππH(ejθ)X(ejθ)dθ若{\quad}y(n)=h(n)x(n)\\ 则有{\quad}Y(e^{j\omega})=\frac{1}{2\pi}H(e^{j\omega})*X(e^{j\omega})\\ =\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\theta})X(e^{j\theta})d_{\theta}

  • 帕斯维尔定理:信号时域的能量和频域的能量的关系。

周期序列的离散傅里叶级数(DFS)

X(k)~=DFS[x(n)~]=n=0N1x(n)~ej2πNkn(1)x(n)~=IDFS[X(k)~]=1Nk=0N1X(k)~ej2πNkn(2)\widetilde{X(k)}=DFS[\widetilde{x(n)}]=\sum_{n=0}^{N-1}\widetilde{x(n)}e^{-j{\frac{2\pi}{N}}kn}{\qquad}(1)\\ \widetilde{x(n)}=IDFS[\widetilde{X(k)}]=\frac{1}{N}\sum_{k=0}^{N-1}\widetilde{X(k)}e^{j{\frac{2\pi}{N}}kn}{\qquad}(2)

基波分量的频率为2π/N{2\pi}/N,幅度是(1/N)X~(1)(1/N){\widetilde{X}}(1),一个周期序列可以用其DFSDFS系数X~(k)\widetilde{X}(k)表示它的频谱分布规律。

周期序列的傅里叶变换表达式(FT)

复指数序列ejω0ne^{j{\omega_0}n}FTFT

首先推导复指数序列ejω0ne^{j{\omega_0}n}的傅里叶变换如下:

X(ejω)=F[ejω0n]=r=+2πδ(ω(ω0+2πr))rX(e^{j\omega})=\mathscr{F}[e^{j{\omega_0}n}]=\sum_{r=-\infin}^{+\infin}2{\pi}\delta\Big({\omega}-({\omega_0}+2{\pi}r)\Big){\quad}r为整数

上式表明:复指数序列ejω0ne^{j{\omega_0}n}FTFT是在ω0+2πr\omega_0+2{\pi}r处的单位冲激函数,强度为2π2\pi

周期序列的FTFT

一般周期序列都能表示成不同复指数序列ejωne^{j{\omega}n}的叠加,因此再根据线性就可以得出该周期序列的傅里叶变换:

x(n)~=a0ejω0n+a1ejω1n+...+amejωmn\widetilde{x(n)}=a_0e^{j{\omega_0}n}+a_1e^{j{\omega_1}n}+...+a_me^{j{\omega_m}n}\\

或者这样考虑,对一般周期序列x(n)~\widetilde{x(n)}展开成DFSDFS,第k次谐波为(X(k)~/N)ej2πNkn(\widetilde{X(k)}/N)e^{j\frac{2\pi}{N}kn},类似于复指数序列的FTFT,其FTFT为:

r=+2πX(k)~Nδ(ω(2πNk+2πr))r\color{blue}\sum_{r=-\infin}^{+\infin}2{\pi}\frac{\widetilde{X(k)}}{N}\delta\Big({\omega}-({\frac{2\pi}{N}k}+2{\pi}r)\Big){\quad}r为整数

上面是第k次谐波的FTFT,那么整个周期序列的FTFT为:

X(ejω)=F[x(n)~]=k=0N1r=+2πX(k)~Nδ(ω(2πNk+2πr))rX(e^{j\omega})={\mathscr{F}}[\widetilde{x(n)}]=\sum_{k=0}^{N-1}{\color{blue}\sum_{r={-\infin}}^{+\infin}2{\pi}\frac{{\widetilde{X(k)}}}{N}\delta\Big({\omega}-({\frac{2\pi}{N}k}+2{\pi}r)\Big){\quad}r为整数}\\

式中,K=0,1,2,...,N1K=0,1,2,...,N-1。如果让k在-\infin++\infin区间变换,上式可简化为:

X(ejω)=k=2πX(k)~Nδ(ω2πNk)X(k)~=n=0N1x(n)~ej2πNknX(e^{j\omega})=\sum_{k={-\infin}}^{\infin}{2\pi}\frac{\widetilde{X(k)}}{N}\delta(\omega-\frac{2\pi}{N}k)\\ 式中{\qquad}{\widetilde{X(k)}}=\sum_{n=0}^{N-1}\widetilde{x(n)}e^{-j{\frac{2\pi}{N}}kn}

对于同一个周期信号,其DTS和FT的模的形状是一样的,不同的是FT用单位冲激函数表示(用带箭头的竖线表示)。因此周期序列的频谱分布用其DFS或者FT表示都可以,但画图时应注意单位冲激函数的画法。

周期序列的FTFT与周期信号FTFT的联系

还记得周期信号fT(t)f_T(t)的傅里叶变换(FT(jω))(F_T(j\omega))是什么吗?(省略掉推导):

FT(jw)=n=2πFnδ(ωnω0)F(ejω)=k=0N1r=+2πF(k)~Nδ(ω(2πNk+2πr))rF_T(jw)=\sum_{n=-\infin}^{\infin}2{\pi}F_n{\delta}(\omega-n\omega_0)\\ F(e^{j\omega})=\sum_{k=0}^{N-1}{\color{blue}\sum_{r={-\infin}}^{+\infin}2{\pi}\frac{{\widetilde{F(k)}}}{N}\delta\Big({\omega}-({\frac{2\pi}{N}k}+2{\pi}r)\Big){\quad}r为整数}\\

上式表明:周期信号的傅里叶变换由无穷多个出现在谐波频率nω0n\omega_0上的冲激函数组成,每一冲激的强度为傅里叶系数FnF_n乘上2π2\pi

看出来周期序列周期信号的傅里叶变换有什么关系了吗?其实两者是一样的:

  • 都是由不同位置的冲激叠加而成;
  • 冲激的强度都是傅氏系数乘上2π2\pi(复指数序列的系数为1)。

但也有不同点,就是冲激位置的规律有点不一样:

  • 周期信号的冲激在谐波频率nω0n\omega_0
  • 周期序列的冲激在2πNk+2πr,r\frac{2{\pi}}{N}k+2{\pi}r,r为整数​处,具有周期性为2π2\pi

为什么会出现周期性呢?并且周期为2π2\pi,请看数字频率与模拟频率的关系与特点

时域离散信号的FT与模拟信号FT之间的关系

X(ejΩT)=1Tk=Xa(jΩjkΩs)X(ejω)=1Tk=Xa(jω2πkT)Ωs=2πFs=2πTX(e^{j{\Omega}T})=\frac{1}{T}\sum_{k=-\infin}^{\infin}X_a(j\Omega-jk\Omega_s)\\ X(e^{j\omega})=\frac{1}{T}\sum_{k=-\infin}^{\infin}X_a(j\frac{\omega-2{\pi}k}{T})\\ \Omega_s=2{\pi}F_s=\frac{2\pi}{T}

上式表明:时域离散信号的频谱也是(前面说的是理想采样信号)模拟信号频谱的周期性延拓,周期为Ωs=2πFs=2πT\Omega_s=2{\pi}F_s=\frac{2{\pi}}{T}

序列的ZZ变换

Z变换的定义及与序列的FT的关系

X(z)=def=n=+x(n)znn=+x(n)zn<X(z)=^{def}=\sum_{n=-\infin}^{+\infin}x(n)z^{-n}\\ \sum_{n=-\infin}^{+\infin}\Big|x(n)z^{-n}\Big|<{\infin}\\

ZZ变换的收敛域一般是一个环带状,即Rx<z<Rx+R_{x-}<\big|z\big|<R_{x+}。如果ZZ变换的收敛域包括单位圆,那么:

X(ejω)=X(z)z=ejωX(e^{j\omega})=X(z)\Big|_{z=e^{j\omega}}

即单位圆上的Z变换就是序列的傅里叶变换(前提是在单位圆上收敛)。

Z变换的性质

X(z)=[x(n)]Rx<z<Rx+X(z)=\big[x(n)\big]{\quad}R_{x-}<\big|z\big|<R_{x+},则有如下性质:

  • 线性

  • 移位性质:

    [x(nn0)]=zn0X(z)Rx<z<Rx+ \big[x(n-n_0)\big]=z^{-n_0}X(z){\quad}R_{x-}<\big|z\big|<R_{x+}

  • 序列乘以指数序列的性质:

    y(n)=anx(n)Y(z)=[anx(n)]=X(a1Z)aRx<z<aRx+ y(n)=a^{n}x(n)\\ Y(z)=\big[a^{n}x(n)\big]=X(a^{-1}Z){\quad}|a|R_{x-}<\big|z\big|<|a|R_{x+}

  • 序列乘以n的ZT

    X(z)=[nx(n)]=zdX(z)dzRx<z<Rx+ X(z)=\big[nx(n)\big]=-z\frac{dX(z)}{dz}{\quad}R_{x-}<\big|z\big|<R_{x+}

  • 复共轭序列的ZT

    X(z)=[x(n)]=X(z)Rx<z<Rx+ X(z)=\big[x^{*}(n)\big]=X^{*}(z^{*}){\quad}R_{x-}<\big|z\big|<R_{x+}

  • 初值定理:x(n)x(n)是因果序列,则:

    x(0)=limzX(z) x(0)=\lim_{z{\to}\infin}X(z)

  • 终值定理:x(n)x(n)是因果序列,其ZZ变换的极点除可以有一个一阶极点在z=1z=1上,其它极点都在单位圆内,则:

    limnx(n)=limz1(z1)X(z) \lim_{n{\to}\infin}x(n)=\lim_{z{\to}1}(z-1)X(z)

  • 时域卷积定理:W(z)W(z)的收敛域就是X(z)X(z)Y(z)Y(z)的公共收敛域。

    W(z)=[x(n)y(n)]=X(z)Y(z)Rw<z<Rw+Rw+=min[Rx+,Ry+]Rw=max[Rx,Ry] W(z)=\big[x(n)^{*}y(n)\big]=X(z)Y(z){\quad}R_{w-}<\big|z\big|<R_{w+}\\ R_{w+}=\min{\big[R_{x+},R_{y+}\big]}\\ R_{w-}=\max{\big[R_{x-},R_{y-}\big]}\\

  • 复卷积定理:不常用

  • 帕斯维尔定理:不常用

利用Z变换分析信号和系统的频响特性

频率响应函数与系统函数

设系统初始状态为0,系统对δ(n)\delta(n)的响应h(n)h(n)的傅里叶变换H(ejω)H(e^{j\omega})称为系统的频率响应函数,即:

H(ejω)=n=+h(n)ejωn=H(ejωn)ejϕ(ω){H(ejωn)ϕ(ω)H(e^{j\omega})=\sum_{n=-\infin}^{+\infin}h(n)e^{-j{\omega}n}=\Big|H(e^{j{\omega}n})\Big|e^{j\phi(\omega)}\\ \begin{cases} 幅频特性:\Big|H(e^{j{\omega}n})\Big|\\ 相频特性:\phi(\omega) \end{cases}

将序列h(n)h(n)进行Z变换得到H(z)H(z),这就是该系统的系统函数,即:

H(z)=Y(z)X(z)=i=0Mbizii=0NaiziH(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{i=0}^{M}b_iz^{-i}}{\sum_{i=0}^{N}a_iz^{-i}}

如果Z变换的收敛域包含单位圆,则:

H(ejω)=H(z)z=ejωH(e^{j\omega})=H(z)\Big|_{z=e^{j\omega}}

输入序列ejωne^{j{\omega}n}的频率响应

y(n)=h(n)x(n)=...=H(ejω)ejωn=H(ejωn)ej[ωn+ϕ(ω)]y(n)=h(n)^{*}x(n)=...=H(e^{j\omega})e^{j{\omega}n}\\ =\Big|H(e^{j{\omega}n})\Big|e^{j\big[{\omega}n+\phi(\omega)\big]}\\

上式说明:单频复指数序列ejωne^{j{\omega}n}通过频率响应函数为H(ejω)H(e^{j\omega})的系统后,输出还是单频复指数序列,只不过幅度放大H(ejωn)\Big|H(e^{j{\omega}n})\Big|倍,相移为ϕ(ω)\phi(\omega)

利用系统的零极点分布分析系统的频率响应特性

H(z)=Ar=1M(1crz1)r=1N(1drz1)H(ejω)=Ar=1Mzrr=1NPrϕ(ω)=ω(NM)+r=1Nαrr=1MβrH(z)=A\frac{\prod_{r=1}^M(1-c_rz^{-1})}{\prod_{r=1}^{N}(1-d_rz^{-1})} \\ |H(e^{j\omega})|=|A|\frac{\prod_{r=1}^M|z_r|}{\prod_{r=1}^N{|P_r|}} \\ \phi(\omega)=\omega(N-M)+\sum_{r=1}^{N}\alpha_r-\sum_{r=1}^{M}\beta_r

即:

  • 幅频响应等于所有的零点到ω\omega的模长的乘积除以所有的极点到ω\omega的模长的乘积;
  • 相频响应等于所有的零点与ω\omega形成的夹角之和减去所有的极点与ω\omega形成的夹角之和。

几种特殊系统的系统函数及特点

全通滤波器

H(ejω)=10ω2π|H(e^{j\omega})|=1{\quad}0{\le}\omega{\le}2\pi

梳妆滤波器

H(zN)=1zN1azN0<a<1H(z^{N})=\frac{1-z^{-N}}{1-az^{-N}}{\qquad}0<a<1

梳妆滤波器可以滤除输入信号中ω=2πNk, k=0,1,...,N1\omega=\frac{2\pi}{N}k,\ k=0,1,...,N-1的频率分量,可以滤除电网谐波干扰和其它频谱等间隔分布的干扰。当a=1a=1时就变成了全通滤波器,下面给出频率响应和零极点分布的MATLABMATLAB代码:

% y(n) -0.9y(n-8) = x(n) - x(n-8);
% H(z) = (1 - z^-8) / (1 - 0.9z^-8)
B = [1 0 0 0 0 0 0 0 -1];
A = [1 0 0 0 0 0 0 0 -0.9];
% impz(B, A, 200);
% [Z, P, K] = tf2zp(B, A)
zplane(B, A);
legend('零点', '极点');
title('$$ H(z) = \frac{1-z^{-8}}{1-0.9z^{-8}}$$', ...
    'Interpreter', 'latex');

[H, w] = freqz(B, A, 500, 'whole');
Hm = abs(H);
Hp = angle(H);
subplot(2, 1, 1);
plot(w, Hm), grid on;
xlabel('\omega(rad/s)');
ylabel('Magnitude');
title('幅频特性曲线');
subplot(2, 1, 2);
plot(w, Hp), grid on;
xlabel('\omega(rad/s)');
ylabel('Phase');
title('相频特性曲线');

a=0.2a=0.2时的频率响应和零极点分布图:

a=0.9a=0.9时的频率响应和零极点分布图:

离散傅里叶变换(DFT)

对离散时间序列x(n)x(n)进行FT得到的结果在频域仍然是连续的,不便于计算机处理。因此目的是对频域也进行离散化,方法有多种,如DFT、DFS。对于DFS,主要步骤就是对x(n)x(n)进行周期性延拓,然后计算其傅里叶级数,该结果和DFT的结果是有关联的。

DFT的定义

X(k)=DFT[x(n)]=n=0N1x(n)WNknk=0,1,...,N1x(n)=IDFT[X(k)]=1Nk=0N1X(K)WNknn=0,1,...,N1WN=ej2πNX(k)=DFT[x(n)]=\sum_{n=0}^{N-1}x(n)W_N^{kn}{\quad}k=0,1,...,N-1\\ x(n)=IDFT[X(k)]=\frac{1}{N}\sum_{k=0}^{N-1}X(K)W_N^{-kn}{\quad}n=0,1,...,N-1\\ W_N=e^{-j\frac{2\pi}{N}}\\

WN=ej2πNW_N=e^{-j\frac{2\pi}{N}}带入即得:

X(k)=DFT[x(n)]=n=0N1x(n)ej2πNknk=0,1,...,N1x(n)=IDFT[X(k)]=1Nk=0N1X(k)ej2πNknn=0,1,...,N1X(k)=DFT[x(n)]=\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn}{\quad}k=0,1,...,N-1\\ x(n)=IDFT[X(k)]=\frac{1}{N}\sum_{k=0}^{N-1}X(k)e^{j\frac{2\pi}{N}kn}{\quad}n=0,1,...,N-1

和DFS的比较

DFS:{X(k)~=DFS[x(n)~]=n=0N1x(n)~ej2πNkn(1)x(n)~=IDFS[X(k)~]=1Nk=0N1X(k)~ej2πNkn(2)DFT:{X(k)=DFT[x(n)]=n=0N1x(n)ej2πNknk=0,1,...,N1(3)x(n)=IDFT[X(k)]=1Nk=0N1X(k)ej2πNknn=0,1,...,N1(4)DFS: \begin{cases} \widetilde{X(k)}=DFS[\widetilde{x(n)}]=\sum_{n=0}^{N-1}\widetilde{x(n)}e^{-j{\frac{2\pi}{N}}kn}{\qquad}(1)\\ \\ \widetilde{x(n)}=IDFS[\widetilde{X(k)}]=\frac{1}{N}\sum_{k=0}^{N-1}\widetilde{X(k)}e^{j{\frac{2\pi}{N}}kn}{\qquad}(2)\\ \end{cases} \\ DFT: \begin{cases} X(k)=DFT[x(n)]=\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn}{\quad}k=0,1,...,N-1{\quad}(3)\\ \\ x(n)=IDFT[X(k)]=\frac{1}{N}\sum_{k=0}^{N-1}X(k)e^{j\frac{2\pi}{N}kn}{\quad}n=0,1,...,N-1{\quad}(4)\\ \end{cases}

其实没有什么不同,只是对序列的取值有限定,这样的限定有如下结论:

  • 如果s(n)s(n)是周期序列,周期为NN,那么其一个周期s1N(n)s_{1N}(n)NNDFTDFT就等于其DFSDFS的主值序列。即:

设序列x(n)x(n)的长度为M。

X(z)=ZT[x(n)]=n=0M1x(n)zn(1)X(ejω)=DTFT[x(n)]=n=0M1x(n)ejωn(2)X(k)=DFT[x(n)]N=n=0M1x(n)WNkn=n=0M1x(n)ej2πNnkk=0,1,...,N1(3)X(z)=ZT[x(n)]=\sum_{n=0}^{M-1}x(n)z^{-n}{\quad}(1)\\ X(e^{j\omega})=DTFT[x(n)]=\sum_{n=0}^{M-1}x(n)e^{-j{\omega}n}{\quad}(2)\\ X(k)=DFT[x(n)]_N=\sum_{n=0}^{M-1}x(n)W_N^{kn}\\ =\sum_{n=0}^{M-1}x(n)e^{-j\frac{2\pi}{N}nk}{\quad}k=0,1,...,N-1{\quad}(3)

  • 序列x(n)x(n)的N点DFT-X(k)X(k)x(n)x(n)的Z变换在单位圆上的N点等间隔采样;

WNk=WNk+mNk,mNX(k+mN)=...=X(k){\because}{\quad}W_N^k=W_N^{k+mN}{\qquad}k,m为整数,N为自然数\\ {\therefore}{\quad}X(k+mN)=...=X(k)

任何周期为N的周期序列x(n)~\widetilde{x(n)}都可以看做长度为N的有限长序列x(n)x(n)的周期性延拓序列,而x(n)x(n)x(n)~\widetilde{x(n)}的一个周期,即:

x(n)~=m=+x(n+mN)x(n)=x(n)~RN(n)\widetilde{x(n)}=\sum_{m=-\infin}^{+\infin}x(n+mN)\\ x(n)=\widetilde{x(n)}R_N(n)

有限长序列x(n)x(n)的N点DFT的结果X(k)X(k)恰好是x(n)x(n)的周期性延拓序列x((n))Nx((n))_N的DFS的结果X(k)~\widetilde{X(k)}的主值序列,即:

X(k)=X(k)~RN(k)X(k)=\widetilde{X(k)}R_N(k)

因此X(k)X(k)实质上就是x(n)x(n)的周期延拓序列x((n))Nx((n))_N的频谱特性。

DFT的性质

线性

时/频域循环移位定理

{y(n)=x((n+m))NRN(n)Y(k)=DFT[y(n)]=WNkmX(k){Y(k)=X((k+l))NRN(k)y(n)=IDFT[Y(K)]N=WNnlx(n)时域循环移位定理 \begin{cases} & y(n)=x((n+m))_NR_N(n)\\ & Y(k)=DFT[y(n)]=W_N^{-km}X(k)\\ \end{cases} \\ 频域循环移位定理 \begin{cases} & Y(k)=X((k+l))_NR_N(k)\\ & y(n)=IDFT[Y(K)]_N=W_N^{nl}x(n) \end{cases}

循环卷积定理

ifx(n)=x2(n)Lx1(n)=[m=0N1x2(m)x1((nm))N]RN(n)haveX(k)=DFT[x(n)]N=X1(k)X2(k)X1(k)=DFT[x1(n)]X2(k)=DFT[x2(n)]if{\quad}x(n)=x_2(n){\color{blue}*_L}x_1(n)\\ =[\sum_{m=0}^{N-1}x_2(m)x_1((n-m))_N]R_N(n) \\ have{\quad}X(k)=DFT[x(n)]_N=X_1(k)X_2(k)\\ X_1(k)=DFT[x_1(n)]\\ X_2(k)=DFT[x_2(n)]\\

复共轭序列的DFT及其对称性

x(n)=xep(n)+xop(n)X(k)=XR(k)+jXI(k)XR(k)=DFT[xep(n)]jXI(k)=DFT[xop(n)]x(n)=x_{ep}(n)+x_{op}(n)\\ X(k)=X_R(k)+jX_I(k)\\ X_R(k)=DFT[x_{ep}(n)]\\ jX_I(k)=DFT[x_{op}(n)]\\

线性卷积和循环卷积的关系

yc(n)=[i=yl(n+iL)]RL(n)y_c(n)=\Big[\sum_{i=-\infin}^{\infin}y_l(n+iL)\Big]R_L(n)

循环卷积等于线性卷积以L为周期进行周期延拓后的序列的主值序列,若要求循环卷积和线性卷积相等,则需要周期延拓后无混叠,即LN+M1L{\ge}N+M-1

频率域采样

X(k)X(k)恢复出x(n)x(n)的条件

如何从x(n)N推导出x(n)呢?下面直接给出结论:

xN(n)=x(n)~RN(n)=[i=x(n+iN)]RN(n)x_N(n)=\widetilde{x(n)}R_N(n)\\ =\Big[\sum_{i=-\infin}^{\infin}x(n+iN)\Big]R_N(n)

如果序列x(n)x(n)的长度为M,则只有当频域采样点数NMN{\ge}M时,才有:

xN(n)=IDFT[X(k)]=x(n)x_N(n)=IDFT\Big[X(k)\Big]=x(n)

即可由频域采样值X(k)X(k)恢复出原序列x(n)x(n),否则产生时域混叠现象。

如何由X(k)X(k)恢复X(z)X(z)X(ejω)X(e^{j\omega})

结论如下,省略推导:

X(z)=k=0N1X(k)1N1zN1WNkz1=k=0N1X(k)ϕk(z)X(z)=\sum_{k=0}^{N-1}X(k)\frac{1}{N}\frac{1-z^{-N}}{1-W_N^{-k}z^{-1}}\\ =\sum_{k=0}^{N-1}X(k){\phi}_k(z)\\

z=ejωz=e^{j\omega}带入上式得:

X(ejω)=k=0N1X(k)ϕ(ω2πNk)ϕ(ω)=1Nsin(ωN/2)sin(ω/2)ejω(N12)X(e^{j\omega})=\sum_{k=0}^{N-1}X(k){\phi}({\omega}-\frac{2\pi}{N}k)\\ {\phi}(\omega)=\frac{1}{N}\frac{sin({\omega}N/2)}{\sin({\omega}/2)}e^{-j{\omega}(\frac{N-1}{2})}

上式中X(z)X(z)表示内插公式ϕk(z)\phi_k(z)称为内插函数

用DFT对信号进行频谱分析

对连续信号进行频谱分析

Xa(kF)=TX(k)=TDFT[x(n)N]k=0,1,...,N1X_a(kF)=TX(k)=T{\cdot}DFT[x(n)_N]{\quad}k=0,1,...,N-1

上式表明,可以通过对连续信号采样并进行DFTDFT再乘TT,近似得到模拟信号频谱的周期延拓函数在第一个周期[0,Fs]\Big[0,F_s\Big]上的N点等间隔采样Xa(kF)X_a(kF)

TpT_pNN可以按照下面两式进行选择:

N>2fcFTp>1F N > \frac{2f_c}{F}\\ T_p > \frac{1}{F}

fcf_c是信号最高频率,FF是谱分辨率,TpT_p是信号持续时间。

对序列进行频谱分析

误差分析

混叠现象

采样速率必须满足采样定理,否则会在ω=π\omega=\pi(对应模拟频率f=Fs/2f=F_s/2)附近发生频谱混叠现象。一般在采样前进行预滤波,滤除高于折叠频率Fs/2F_s/2的频率成分,以免发生频谱混叠现象。

栅栏效应

N点DFT是在频率区间[0,2π][0,2\pi]上对时域离散信号的频谱进行N点等间隔采样,而采样点之间的频谱是看不到的。对有限长序列,可以在原序列尾部补零;对无限长序列,可以增大截取长度DFT变换区间长度,从而使频域采样间隔变小,增加频域采样点数和采样点位置,使原来漏掉的某些频谱成分被检测出来。

截断效应

实际中遇到的x(n)x(n)可能是无限长的,用DFT对其进行频谱分析时,必须将其截短,形成有限长序列y(n)=x(n)w(n)y(n)=x(n)w(n)w(n)w(n)称为窗函数,长度为N。

赞赏