精华内容
下载资源
问答
  • 离散时间信号和线性时不变离散时间系统的频域分析

    一、实验目的:

    熟悉Mat lab中的信号处理工具箱及相关命令,掌握用Mat lab验证离散时间傅里叶变换相关性质的方法;掌握用Mat lab分析离散时间系统频域特性的方法并能绘制幅频特性和相频特性图。

    二、实验内容及要求:

    1. 离散时间信号的频域分析:

    (1)修改程序P3.1,在范围 内计算如下序列的离散时间傅里叶变换并绘制幅频特性和相频特性图。
    在这里插入图片描述
    (2)修改程序P3.2,D=8,验证离散时间傅里叶变换的时移特性并绘图。
    (3)修改程序P3.3,w = 0.2π ,验证离散时间傅里叶变换的频移特性并绘图。

    2. 线性时不变离散时间系统的频域分析:

    (1)修改程序P3.1,计算并画出如下传输函数的频率响应,它是哪种类型的滤波器?
    在这里插入图片描述
    (2)修改程序P4.2,计算并画出一个长度为6的滑动平均滤波器的增益响应。
    (3)用Matlab画出如下滤波器的零极点图并分析稳定性。
    在这里插入图片描述

    三、实验结果及问题回答:

    1. 离散时间信号的频域分析:

    (1)实验结果:

    % Program P3_1
    % Evaluation of the DTFT 
    clf;
    % Compute the frequency samples of the DTFT
    w = 0*pi:8*pi/511:pi;
    num = [0.7 -0.5 0.3 1];den = [1 0.3 -0.5 0.7];
    h = freqz(num, den, w);
    % Plot the DTFT
    subplot(2,1,1)
    plot(w/pi,real(h));grid
    title('Real part of H(e^{j\omega})')
    xlabel('\omega /\pi');
    ylabel('Amplitude');
    subplot(2,1,2)
    plot(w/pi,imag(h));grid
    title('Imaginary part of H(e^{j\omega})')
    xlabel('\omega /\pi');
    ylabel('Amplitude');
    pause
    subplot(2,1,1)
    plot(w/pi,abs(h));grid
    title('Magnitude Spectrum |H(e^{j\omega})|')
    xlabel('\omega /\pi');
    ylabel('Amplitude');
    subplot(2,1,2)
    plot(w/pi,angle(h));grid
    title('Phase Spectrum arg[H(e^{j\omega})]')
    xlabel('\omega /\pi');
    ylabel('Phase in radians');
    

    在这里插入图片描述
    (2)实验结果:

    % Program P3_2
    % Time-Shifting Properties of DTFT
    clf;
    w = -pi:2*pi/255:pi; wo = 0.4*pi; D = 8;
    num = [1 2 3 4 5 6 7 8 9];
    h1 = freqz(num, 1, w);
    h2 = freqz([zeros(1,D) num], 1, w);
    subplot(2,2,1)
    plot(w/pi,abs(h1));grid
    title('Magnitude Spectrum of Original Sequence')
    subplot(2,2,2)
    plot(w/pi,abs(h2));grid
    title('Magnitude Spectrum of Time-Shifted Sequence')
    subplot(2,2,3)
    plot(w/pi,angle(h1));grid
    title('Phase Spectrum of Original Sequence')
    subplot(2,2,4)
    plot(w/pi,angle(h2));grid
    title('Phase Spectrum of Time-Shifted Sequence')
    

    在这里插入图片描述
    (3)实验结果:

    % Program P3_3
    % Frequency-Shifting Properties of DTFT
    clf;
    w = -pi:2*pi/255:pi; wo = 0.2*pi;
    num1 = [1 3 5 7 9 11 13 15 17];
    L = length(num1);
    h1 = freqz(num1, 1, w);
    n = 0:L-1;
    num2 = exp(wo*i*n).*num1;
    h2 = freqz(num2, 1, w);
    subplot(2,2,1)
    plot(w/pi,abs(h1));grid
    title('Magnitude Spectrum of Original Sequence')
    subplot(2,2,2)
    plot(w/pi,abs(h2));grid
    title('Magnitude Spectrum of Frequency-Shifted Sequence')
    subplot(2,2,3)
    plot(w/pi,angle(h1));grid
    title('Phase Spectrum of Original Sequence')
    subplot(2,2,4)
    plot(w/pi,angle(h2));grid
    title('Phase Spectrum of Frequency-Shifted Sequence')
    

    在这里插入图片描述

    2. 线性时不变离散时间系统的频域分析:

    (1)实验结果:
    % Program P3_1
    % Evaluation of the DTFT
    clf;
    % Compute the frequency samples of the DTFT
    w = 0pi:8pi/511:1*pi;
    num = [0.15 0 -0.15];den = [1 -0.5 0.7];
    h = freqz(num, den, w);
    % Plot the DTFT
    subplot(2,1,1)
    plot(w/pi,real(h));grid
    title(‘Real part of H(e^{j\omega})’)
    xlabel(’\omega /\pi’);
    ylabel(‘Amplitude’);
    subplot(2,1,2)
    plot(w/pi,imag(h));grid
    title(‘Imaginary part of H(e^{j\omega})’)
    xlabel(’\omega /\pi’);
    ylabel(‘Amplitude’);
    pause
    subplot(2,1,1)
    plot(w/pi,abs(h));grid
    title(‘Magnitude Spectrum |H(e^{j\omega})|’)
    xlabel(’\omega /\pi’);
    ylabel(‘Amplitude’);
    subplot(2,1,2)
    plot(w/pi,angle(h));grid
    title(‘Phase Spectrum arg[H(e^{j\omega})]’)
    xlabel(’\omega /\pi’);
    ylabel(‘Phase in radians’);

    在这里插入图片描述
    (2)实验结果:

    % Program P4_2
    % Gain Response of a Moving Average Lowpass Filter
    clf;
    M = 6;
    num = ones(1,M)/M;
    [g,w] = gain(num,1);
    plot(w/pi,g);grid
    axis([0 1 -50 0.5])
    xlabel('\omega /\pi');ylabel('Gain in dB');
    title(['M = ', num2str(M)])
    
    function [g,w] = gain(num,den)
    % Computes the gain function in dB of a 
    % transfer function at 256 equally spaced points 
    % on the top half of the unit circle
    % Numerator coefficients are in vector num
    % Denominator coefficients are in vector den
    % Frequency values are returned in vector w
    % Gain values are returned in vector g
    w = 0:pi/255:pi;
    h = freqz(num,den,w);
    g = 20*log10(abs(h));
    

    在这里插入图片描述
    (3)实验结果:

    clc;
    clear;
    num = [ 1 0 0];
    den = [1 -1.848 0.85];
    zplane(num,den)

    在这里插入图片描述

    展开全文
  • 在连续时间信号与系统中,信号一般用连续变量时间t函数表示,系统用微分方程描述,其频域分析方法是拉普拉斯变换和傅立叶变换。在时域离散信号与系统中,信号用序列表示,其自变量仅取整数,非整数时无定义,系统...

    信号与系统的分析方法有两种:时域分析方法和频域分析方法。

    在连续时间信号与系统中,信号一般用连续变量时间t的函数表示,系统用微分方程描述,其频域分析方法是拉普拉斯变换和傅立叶变换。在时域离散信号与系统中,信号用序列表示,其自变量仅取整数,非整数时无定义,系统则用差分方程描述,频域分析方法是Z变换和序列傅立叶变换法。

    Z变换在离散时间系统中的作用就如同拉普拉斯变换在连续时间系统中的作用一样,它把描述离散系统的差分方程转化为简单的代数方程,使其求解大大简化。因此,对求解离散时间系统而言,Z变换是一个极重要的数学工具。

    2.2序列的傅立叶变换(离散时间傅立叶变换)

    一、序列傅立叶变换:

    正变换:DTFT[x(n)]=2724041_1.gif(2.2.1)

    反变换:DTFT-12724041_2.gif

    式(2.2.1)级数收敛条件为

    |2724041_3.gif|=2724041_4.gif   (2.2.2)

    上式称为x(n)绝对可和。这也是DTFT存在的充分必要条件。当遇到一些绝对不可和的序列,例如周期序列,其DTFT可用冲激函数的形式表示出来。

    二、序列傅立叶变换的基本性质:

    1、DTFT的周期性

    2724041_1.gif2724041_5.gif是频率w的周期函数,周期为2p。

    2724041_6.gif = 2724041_3.gif

    问题1:设x(n)=RN(n),求x(n)的DTFT。

    2724041_1.gif=2724041_7.gif=2724041_8.gif

    =2724041_9.gif=2724041_10.gif

    设N为4,画出幅度与相位曲线。

    2724041_11.jpg

    2、线性

    2724041_12.gif=DTFT[x1(n)],2724041_13.gif=DTFT[x2(n)],

    则:DTFT[a x1(n)+b x2(n)]

    = 2724041_14.gif = a2724041_12.gif+b2724041_13.gif

    3、序列的移位和频移

    2724041_5.gif = DTFT[x(n)],

    则:DTFT[x(n-n0)] = 2724041_15.gif

    = 2724041_16.gif2724041_5.gif

    DTFT[2724041_17.gifx(n)]  = 2724041_18.gif

    = 2724041_19.gif = 2724041_20.gif

    4、DTFT的对称性

    共轭对称序列的定义:设序列2724041_21.gif满足下式

    2724041_22.gif

    则称2724041_21.gif为共轭对称序列。

    共轭对称序列的性质:

    共轭对称序列的实部是偶函数,虚部是奇函数

    证明:2724041_21.gif=2724041_23.gif+j2724041_24.gif(实部加虚部)

    2724041_22.gif

    2724041_23.gif+j2724041_24.gif=2724041_25.gif-j2724041_26.gif

    2724041_23.gif=2724041_25.gif(偶函数)

    2724041_24.gif=-2724041_26.gif(奇函数)

    一般情况下,共轭对称序列用2724041_27.gif表示:

    2724041_28.gif

    共轭反对称序列的定义:设序列2724041_21.gif满足下式

    2724041_29.gif

    则称2724041_21.gif为共轭反对称序列。

    共轭反对称序列的性质:

    共轭反对称序列的实部是奇函数,虚部是偶函数

    证明:2724041_21.gif=2724041_23.gif+j2724041_24.gif(实部加虚部)

    2724041_29.gif

    2724041_23.gif+j2724041_24.gif=2724041_30.gif+j2724041_26.gif

    2724041_23.gif=2724041_30.gif(奇函数)

    2724041_24.gif=2724041_26.gif(偶函数)

    一般情况下,用2724041_31.gif来表示

    2724041_32.gif

    一个序列可用共轭对称序列2724041_27.gif与共轭反对称序列2724041_31.gif之和表示。即:

    x(n)= 2724041_27.gif+ 2724041_31.gif    (2.2.16)

    问题1:2724041_33.gif =?

    2724041_33.gif=2724041_34.gif+2724041_35.gif

    =2724041_27.gif2724041_31.gif

    2724041_33.gif= 2724041_27.gif2724041_31.gif    (2.2.17)

    2724041_27.gif=2724041_36.gif(2724041_37.gif+2724041_33.gif)

    2724041_31.gif=2724041_36.gif(2724041_37.gif2724041_33.gif)

    对于频域函数2724041_5.gif,也可分解成共轭对称分量和共轭反对称分量之和:

    2724041_38.gif

    式中,2724041_39.gif是共轭对称分量,2724041_40.gif是共轭反对称分量,它们满足:

    2724041_39.gif=2724041_41.gif2724041_40.gif=2724041_42.gif

    且:2724041_43.gif

    2724041_44.gif

    2724041_39.gif:共轭对称分量,它的实部是偶函数,虚部是奇函数;2724041_40.gif:共轭反对称分量,它的实部是奇函数,虚部是偶函数。

    下面研究DTFT的对称性,按下面两部分进行分析

    a)将序列x(n)分成实部与虚部,即:

    2724041_21.gif=2724041_23.gif+j2724041_24.gif(2724041_23.gif2724041_24.gif都是实数序列)

    则:2724041_38.gif

    式中:2724041_39.gif=DTFT[2724041_23.gif]=2724041_45.gif

    2724041_40.gif=DTFT[j2724041_24.gif]=j2724041_46.gif

    结论:序列分成实部与虚部两部分,实部对应于2724041_5.gif中的2724041_39.gif,虚部和j一起对应于2724041_5.gif中的2724041_40.gif

    b)将序列分成共轭对称部分2724041_27.gif和共轭反对称部分2724041_31.gif,x(n)= 2724041_27.gif+ 2724041_31.gif

    2724041_27.gif=2724041_36.gif(2724041_37.gif+2724041_33.gif)

    2724041_31.gif=2724041_36.gif(2724041_37.gif2724041_33.gif)

    将上面两式分别进行DTFT,得到:

    DTFT[2724041_27.gif]=2724041_36.gif(2724041_5.gif+2724041_47.gif)=Re[2724041_5.gif]=2724041_48.gif

    DTFT[2724041_31.gif]=2724041_36.gif(2724041_5.gif2724041_49.gif)=jIm[2724041_5.gif]=j2724041_50.gif

    2724041_5.gif=2724041_48.gif+j2724041_50.gif

    x(n)= 2724041_27.gif+ 2724041_31.gif

    结论:序列的共轭对称部分2724041_27.gif对应于2724041_5.gif的实部2724041_48.gif,而序列的共轭反对称部分2724041_31.gif对应于2724041_5.gif的虚部加j。

    应用:利用DTFT的对称性讨论当h(n)是实序列时,其DTFT的特性。

    ∵h(n)是实序列,所以它所对应的DTFT:2724041_51.gif=2724041_52.gif,具有共轭对称性,2724041_51.gif的实部偶对称,虚部奇对称。

    5、时域卷积定理

    设y(n)=x(n)*h(n)

    则:2724041_53.gif=2724041_5.gif×2724041_51.gif=2724041_5.gif2724041_51.gif  (2.2.32)

    证明:y(n)= x(n)*h(n)=2724041_54.gif

    2724041_53.gif=DTFT[y(n)]

    = 2724041_55.gif=2724041_56.gif

    =2724041_57.gif

    =2724041_58.gif

    =2724041_51.gif2724041_5.gif

    6、频域卷积定理

    设y(n) = x(n) h(n)

    2724041_53.gif=2724041_59.gif2724041_5.gif*2724041_51.gif=2724041_59.gif2724041_60.gif

    =2724041_59.gif2724041_61.gif

    证明:2724041_53.gif=2724041_55.gif=2724041_62.gif

    =2724041_63.gif

    =2724041_59.gif2724041_64.gif2724041_65.gif

    =2724041_59.gif2724041_64.gif2724041_66.gif

    =2724041_59.gif2724041_5.gif*2724041_51.gif

    7、Parseval(帕斯维尔)(帕塞瓦尔)定理

    2724041_67.gif=2724041_59.gif2724041_68.gif (2.2.34)

    证明:

    2724041_67.gif=2724041_69.gif=2724041_70.gif

    = 2724041_59.gif2724041_71.gif2724041_72.gif

    =2724041_59.gif2724041_73.gif

    =2724041_59.gif2724041_68.gif

    2.5  Z变换的定义与收敛域

    一、Z变换的定义

    若序列为x(n),则幂级数

    2724041_74.gif(2.5.1)

    称为序列x(n)的Z变换,也称为双边Z变换。式中z为复变量,它所在的复平面称为z平面。亦可将x(n)的Z变换表示为

    ZT[x(n)] = X(z)

    二、Z变换的收敛域

    我们知道,2724041_74.gif是一幂级数,只有收敛时Z变换才有意义。X(z)收敛的条件是:

    2724041_75.gif      (2.5.3)

    X(z)能够收敛的z取值集合称为X(z)的收敛域。

    一般收敛域用环状域表示。即:

    2724041_76.gif

    2724041_77.jpg

    ∴Z变换的公式

    2724041_74.gif   2724041_76.gif(2.5.1)

    常见的Z变换是一个有理函数,表示为:

    2724041_78.gif

    分子多项式2724041_79.gif的根是2724041_80.gif的零点,分母多项式2724041_81.gif的根是2724041_80.gif的极点。在极点处Z变换不存在。因此收敛域中没有极点,收敛域总是用极点限定其边界。

    1、有限长序列Z变换的收敛域

    有限长序列是指在有限区间n1≤n≤n2之间序列具有非零的有限值,在此区间外,序列值皆为零。有限长序列Z变换为2724041_82.gif,所以收敛域为

    0

    如n1≥0,收敛域为0

    如n2≤0,收敛域为0≤|z|

    2、右边序列Z变换的收敛域

    右边序列是指在n≥n1时,x(n)有值,在n<n1时,x(n)=0。其Z变换为

    2724041_83.gif

    2724041_84.gif

    此式右端第一项为有限长序列的Z变换,它的收敛域为0≤|z|。综合此两项,只有两项都收敛时级数才收敛。所以右边序列Z变换的收敛域为2724041_86.gif

    因果序列是最重要的一种右边序列,即n1=0的右边序列。收敛域为2724041_85.gif(也可以写成2724041_87.gif),所以,|z|=∞处Z变换收敛是因果序列的特征。

    3、左边序列Z变换的收敛域

    左边序列是指在n≤n2时,x(n)有值,n>n2时,x(n)=0。其Z变换为2724041_88.gif

    2724041_89.gif

    此式第二项是有限长序列的Z变换,收敛域为0

    4、双边序列Z变换的收敛域

    这类序列是指n为任意值时x(n)皆有值的序列。2724041_90.gif双边序列的收敛域为2724041_76.gif

    问题1:求序列x(n)= RN(n)的Z变换及收敛域,并画出收敛域。

    解:X(z)=2724041_91.gif=2724041_92.gif。因为这是有限长序列,所以收敛域为0

    思考:RN(n)的DTFT存在吗?

    问题2:x(n)=anu(n),求其Z变换及收敛域,并画出收敛域。

    解:这是右边序列,且是因果序列,其Z变换为X(z)=2724041_93.gif。收敛域为2724041_94.gif(或写成2724041_95.gif)

    思考:anu(n)的DTFT存在吗?

    问题3:x(n)=-anu(-n-1),求其Z变换及收敛域,并画出收敛域。

    解:这是一个左边序列。其Z变换为

    2724041_96.gif

    2724041_97.gif

    收敛域为0≤|z|

    思考:-anu(-n-1)的DTFT存在吗?

    结论:当Z变换的收敛域中包含单位圆时,用Z变换可求出DTFT。

    2724041_5.gif=2724041_98.gif(2.5.4)

    上式称为单位圆上的Z变换就是离散时间傅立叶变换。

    回顾:观察零极点。

    结论:零点可以在复平面的任意处,但极点在收敛域的边缘或收敛域的外面。

    2.5.3  Z反变换

    已知序列的Z变换及其收敛域,求序列称为Z反变换。表示为x(n)=ZT-1[X(z)]

    2724041_74.gif         2724041_76.gif

    2724041_99.jpg

    2724041_100.gif

    其中,c是X(z)收敛域中一条逆时针的闭合曲线。

    求Z反变换的方法通常有三种:围线积分法(留数法)、部分分式展开法和长除法。

    一、围线积分法(留数法)

    直接计算围线积分比较麻烦,一般都采用留数定理来求解。按留数定理,若函数F(z)=X(z)zn-1在围线c上连续,在c以内有K个极点zk,则有

    2724041_101.jpg(2.5.6)

    设zr是X(z)zn-1的单极点,则根据留数定理:

    2724041_102.gif

    如果zk是L阶极点,则根据留数定理,

    2724041_103.gif

    (2.5.8)

    (2.5.8)表明,对于L阶极点,需要求L-1次导数,这是比较麻烦的。如果c内有多阶极点,而c外没有多阶极点时,可根据留数辅助定理改求c外所有极点之和,使问题简单化。

    若函数F(z)=X(z)zn-1在围线c上连续,在c以内有K个极点zk,而在c以外有M个极点zm,(K,M为有限值)。现在c内有多阶极点,而c外没有多阶极点,根据留数辅助定理改求c外所有极点之和。得:

    2724041_104.jpg

    (2.5.9)(2.5.9)应用条件是X(z)zn-1在z=∞有两阶或二阶以上零点,即要分母多项式z的阶次比分子多项式z的阶次高二阶或二阶以上。

    问题1:已知X(z)=z2/[(4-z)(z-1/4)],1/4

    2724041_105.jpg求Z反变换。

    解:c,c为X(z)的收敛域

    内的闭合围线,画出收敛域及c。

    X(z)zn-1=2724041_106.gif。现在来看极点在围线c内部及外部的分布情况及极点阶数。

    2724041_107.gif时,函数在围线c内只有z=1/4处一个一阶极点,

    2724041_108.gif

    =2724041_109.gif2724041_107.gif

    2724041_110.gif时,函数2724041_106.gif在围线外部只有一个一阶极点z=4,而在围线的内部则有z=1/4处一阶极点及z=0处一(n+1)阶极点,所以采用围线外部的极点较方便。

    2724041_111.gif

    =2724041_112.gif2724041_110.gif

    2724041_113.gif

    问题2:已知X(z)=z2/[(4-z)(z-1/4)],|z|>4,

    求Z反变换。

    2724041_105.jpg解:c,c为X(z)的收敛域

    内的闭合围线。

    X(z)zn-1=2724041_106.gif。现在来看在围c内部及外部的分布情况及极点阶数。

    2724041_107.gif时,

    函数2724041_106.gif在围线c内z=1/4处有一个一阶极点,z=4处有一个一阶极点,

    2724041_114.gif+2724041_115.gif

    =2724041_116.gif2724041_107.gif

    当n=-1时,x(n)=0,∴x(n)= 2724041_116.gif2724041_117.gif

    2724041_110.gif时,函数2724041_106.gif在围线外部没有一个极点,所以采用围线外部的极点较方便。由于围线外部没有一个极点,∴x(n)=0。

    ∴x(n)=(2724041_116.gif)u(n)

    二、部分分式展开法

    对于大多数单极点的序列,常常用这种部分分式展开法求Z反变换。

    X(z)=B(z)/A(z)= X1(z)+ X2(z)+…+ XK(z),则

    2724041_37.gif=ZT-1[X1(z)]+ZT-1[ X2(z)]+…+ZT-1[XK(z)]

    ZT-1[X1(z)]、ZT-1[ X2(z)]、…ZT-1[XK(z)]可从Z变换表中直接查表得出

    问题1:设X(z)=z2/[(z-2)(z-0.5)],|z|>2,

    求Z反变换。

    解:X(z) =z2/[(z-2)(z-0.5)]

    2724041_118.gif

    A1=2724041_119.gif,A2=2724041_120.gif

    2724041_121.gif2724041_122.gif

    ∵收敛域为|z|>2,∴x(n)=2724041_123.gif

    三、幂级数展开法

    因为2724041_37.gif的Z变换定义为z-1的幂级数,即

    2724041_124.gif

    所以只要在给定得收敛域内,把X(z)展成幂级数,则级数的系数就是序列2724041_37.gif

    当X(z)的收敛域为|z|>Rx-时,则2724041_37.gif必为因果序列,此时应将X(z)展成z的负幂级数,为此,X(z)的分子分母应按z的降幂排列;

    当X(z)的收敛域为|z|必为左边序列,X(z)的分子分母应按z的升幂排列;

    问题1:已知2724041_125.gif,|z|>3

    解:因为收敛域|z|>3,所以这是因果序列,因此,X(z)分子分母按z的降幂排列。

    2724041_126.gif

    进行长除

    2.5.4  Z变换的基本性质和定理

    一、线性

    线性就是要满足比例性和可加性。若

    X(z) =ZT [x(n) ],2724041_76.gif

    Y(z) = ZT [y(n) ],2724041_127.gif

    则ZT[ax(n)+by(n)]=aX(z)+b Y(z),2724041_128.gif

    2724041_129.gif2724041_130.gif

    二、序列的移位

    若X(z) = ZT [x(n) ],2724041_76.gif

    则有ZT [x(n-m) ] =z-mX(z),2724041_76.gif

    三、乘以指数序列

    若X(z) =ZT [x(n) ],2724041_76.gif

    则ZT [anx(n) ]=X(2724041_131.gif),2724041_132.gif

    四、序列乘以n

    若X(z) =ZT [x(n) ],2724041_76.gif

    则ZT [n x(n) ]=-z2724041_133.gif2724041_76.gif

    五、复序列取共扼

    一个复序列x(n)的共扼序列为x*(n)

    若ZT [x(n) ] =X(z),2724041_76.gif

    则ZT [x*(n) ] =X*(z*),2724041_76.gif

    六、翻转序列

    若ZT [x(n) ] =X(z),2724041_76.gif

    则ZT [x(-n) ] =X(2724041_134.gif),2724041_135.gif

    七、(因果序列)初值定理

    对于因果序列x(n),即x(n)=0,n<0,ZT[x(n) ] =X(z)有

    2724041_136.gif

    八、(因果序列)终值定理

    设x(n)为因果序列,且X(z) = ZT [x(n) ]的极点处于单位圆|z|=1以内(单位圆上最多在z=1处可有一阶极点),则

    2724041_137.gif

    九、序列的卷积和(时域卷积和定理)

    设y(n)为x(n)与h(n)的卷积和

    y(n)= x(n)*h(n)=2724041_54.gif

    X(z) =ZT [x(n) ],2724041_76.gif

    H(z) =ZT [h(n) ],2724041_138.gif

    则Y(z) =ZT [y(n) ]= X(z) H(z),

    2724041_139.gif

    十、序列相乘(z域卷积定理)

    若y(n)= x(n)·h(n),且

    X(z) =ZT [x(n) ],2724041_76.gif

    H(z) =ZT [h(n) ],2724041_138.gif

    则Y(z) =ZT [y(n) ]=ZT[x(n)·h(n)]

    =2724041_140.gif2724041_141.gif

    其中c是v平面上,2724041_142.gif与H(v)的公共收敛域内环绕原点的一条反时针旋转的单封闭围线。

    v平面收敛域为

    2724041_143.gif

    或Y(z) =ZT [y(n) ]=ZT[x(n)·h(n)]

    =2724041_144.gif2724041_141.gif

    其中c是v平面上,2724041_145.gif与X(v)的公共收敛域内环绕原点的一条反时针旋转的单封闭围线。

    v平面收敛域为

    2724041_146.gif

    十一、帕斯维尔(Parseval)定理

    若X(z) =ZT [x(n) ],2724041_76.gif

    H(z) =ZT [h(n) ],2724041_138.gif

    2724041_147.gif

    2724041_148.gif

    v平面上,c所在的收敛域为

    2724041_149.gif

    证明:Y(z) =ZT[x(n)·h*(n)]

    =2724041_150.gif

    =2724041_151.gif2724041_141.gif

    因为2724041_147.gif,所以z=1在收敛域中。令z=1代入上式,

    2724041_152.gif=2724041_148.gif

    v平面上,c所在的收敛域为

    2724041_149.gif

    如果X(z),H(z)在单位圆上都收敛,则c可取为单位圆,即2724041_153.gif,则2724041_154.gif

    如果h(n)=x(n),则进一步有2724041_155.gif

    2.5.5利用Z变换解差分方程

    在第一章中介绍了差分方程的递推解法,下面介绍Z变换解法。这种方法将差分方程变成了代数方程,使求解过程简单

    设N阶线性常系数差分方程为

    2724041_156.gif     (2.5.30)

    一、求2724041_157.gif2724041_158.gif

    对(2.5.30)求双边Z变换:

    2724041_159.gif=2724041_160.gif

    2724041_157.gif=2724041_158.gif/2724041_80.gif=2724041_161.gif,h(n)=ZT-1[2724041_157.gif]

    2724041_158.gif=2724041_157.gif2724041_80.gif,y(n)=ZT-1[2724041_158.gif]

    2.6离散系统的系统函数,系统的频率响应

    信号和系统的频率特性一般用序列的傅立叶变换和Z变换进行分析。

    一、传输函数与系统函数

    设系统初始状态为零,输出端对输入为单位抽样序列d(n)的响应,称为系统的单位抽样响应h(n)。对h(n)进行傅立叶变换得到:2724041_162.gif=2724041_163.gif,一般称为2724041_162.gif为系统的传输函数,它表征系统的频率特性。

    将h(n)进行Z变换,得到2724041_164.gif,一般称H(z)为系统的系统函数,它表征了系统的复频域特性。

    如已知系统的N阶线性常系数差分方程,进行双边Z变换,得到系统函数的一般表示式:

    2724041_157.gif2724041_165.gif

    如果2724041_157.gif的收敛域包含单位圆|z|=1则,2724041_162.gif2724041_157.gif的关系:2724041_162.gif=2724041_166.gif。即单位圆上的系统函数就是系统的传输函数2724041_162.gif

    二、用系统函数的极点分布分析系统的因果性和稳定性

    因果(可实现)系统其单位脉冲响应h(n)一定满足:当n<0时,h(n)=0,那么其系统函数2724041_157.gif的收敛域一定包含¥点。

    系统稳定要求2724041_167.gif,对照ZT定义,系统稳定要求收敛域包含单位圆。

    所以系统因果且稳定,收敛域包含¥点和单位圆,那么收敛域表示为:r

    问题1:一个因果系统的系统函数为2724041_157.gif=2724041_168.gif,其中a为实数,问:a在哪些范围内才能使系统稳定?

    解:因为系统因果,所以收敛域为|a|

    三、利用系统的零极点分布分析系统的频率特性

    2724041_157.gif=2724041_161.gif

    将上式因式分解,得到:

    2724041_157.gif=A2724041_169.gif

    式中,2724041_170.gif2724041_157.gif的零点,2724041_171.gif是其极点。A参数影响传输函数的幅度大小,影响系统特性的是零点和极点的分布。下面采用几何方法研究系统零极点分布对系统频率特性的影响。

    2724041_157.gif=A2724041_172.gif

    设系统稳定,将z=2724041_173.gif代入,得:

    2724041_162.gif=A2724041_174.gif

    在z平面上,2724041_173.gif-2724041_170.gif用一根由零点2724041_170.gif指向单位圆上2724041_173.gif点B的向量2724041_175.gif表示。同样2724041_173.gif-2724041_171.gif用由极点指向2724041_173.gif点B的向量2724041_176.gif表示,如图2.6.2。

    2724041_177.jpg

    将向量用极坐标表示:2724041_175.gif=2724041_170.gif2724041_178.gif2724041_176.gif=2724041_171.gif2724041_179.gif,得到:

    2724041_180.gif=A2724041_181.gif=2724041_182.gif

    2724041_183.gif= |A|2724041_184.gif   (2.6.8)

    2724041_185.gif=2724041_186.gif2724041_187.gif(N=M)(2.6.9)

    当频率w从零变化到2p时,这些向量的终点B沿单位圆逆时针旋转一周,按照(2.6.8)和(2.6.9)分别估算出系统的幅度特性和相位特性。

    按照(2.6.8)知道零极点的分布后,可以很容易地确定零极点位置对系统特性的影响。当B点转到极点附近时,极点矢量长度最短因而幅度特性可能出现峰值,且极点愈靠近单位圆,极点矢量长度愈短,峰值愈高愈尖锐。如果极点在单位圆上,则幅度特性为¥,系统不稳定。对于零点,情况相反,当B点转到零点附近,零点矢量长度变短,幅度特性将出现谷值,零点愈靠近单位圆,谷值愈接近零。当零点处在单位圆上时,谷值为零。总结以上结论:极点位置主要影响频响的峰值位置及尖锐程度,零点位置主要影响频响的谷点位置及形状。1.5模拟信号数字处理方法

    2724041_188.jpg

    数字信号处理技术优于模拟信号处理技术,故人们将模拟信号数字化,即经过采样、量化编码最终形成数字信号。

    连续时间信号变为离散时间信号是由“采样”这一过程完成的。采样是将模拟信号数字化的第一个环节。它是利用周期性抽样脉冲序列(常用p(t)表示)从连续信号中抽取一系列的离散值来得到抽样信号的。如下图,根据每个脉冲宽度的不同,可将抽样分为两种:

    理想抽样实际抽样2724041_189.jpg

    我们要研究的是,信号被抽样后其频谱将会有什么变化?在什么条件下,可从抽样信号2724041_190.gif中不失真地恢复原来信号xa(t)?

    设:xa(t)的傅立叶变换为:2724041_191.gif,抽样脉冲序列p(t)的傅立叶变换为:2724041_192.gif,抽样信号2724041_190.gif的傅立叶变换为:2724041_193.gif

    2724041_190.gif= xa(t)p(t),∴2724041_194.gif

    2724041_195.gif2724041_196.gif

    由上图得,抽样脉冲序列p(t)的周期为T,则抽样频率2724041_197.gif。则周期信号p(t)的傅立叶变换2724041_198.gif,其中2724041_199.gif

    2724041_200.gif

    2724041_190.gif的傅立叶变换为:2724041_193.gif

    一、理想抽样

    p(t)=2724041_201.gif

    2724041_202.gif=2724041_203.gif

    2724041_198.gif=2724041_204.gif

    2724041_200.gif

    =2724041_205.gif

    2724041_190.gif= xa(t) p(t),p(t)=2724041_206.gif

    2724041_190.gif= xa(t) p(t)= xa(t) 2724041_206.gif

    =2724041_207.gif=2724041_208.gif

    2724041_209.gif

    2724041_210.jpg

    2724041_211.jpg

    2724041_197.gif

    二、抽样定理

    要想抽样后能不失真的还原出原信号,则抽样频率必须大于两倍的信号谱最高频率即2724041_212.gif,这就是抽样定理。

    对连续信号进行等间隔采样形成采样信号,采样信号的频谱是原连续信号的频谱以采样频率为周期进行周期性的延拓形成的。

    2724041_213.gif = 2724041_205.gif   (1.5.5)

    三、抽样的恢复

    如果满足抽样定理,则抽样后不会产生频谱混叠,故将2724041_213.gif通过如图所示的理想低通滤波器,就可得到信号频谱。

    2724041_214.jpg

    虽然理想低通滤波器是不可实现的,但在一定精度范围内可以用可实现的滤波器来逼近

    下面讨论如何由抽样值来恢复原来的模拟信号。即2724041_190.gif通过H(jW)系统的响应特性。理想低通滤波器的冲激响应为

    2724041_215.gif2724041_190.gif与h(t)的卷积积分,即得理想低通滤波器的输出为2724041_216.gif

    2724041_217.gif

    2724041_218.gif

    这就是内插值公式,即由信号的抽样值2724041_219.gif经此公式而得到连续信号2724041_220.gif,而2724041_221.gif称为内插函数,如图所示,在抽样点mT上,函数值为1,其余抽样点上,函数值为0。

    2724041_222.jpg

    2724041_223.jpg在每个抽样点上,只有该点所对应的内插函数不为零,这使得各抽样点上信号值不变,而抽样点之间的信号则由各加权抽样函数波形的延伸叠加而成,如下图所示。这个公式说明了只要抽样频率高于两倍信号最高频率,则整个连续信号就可完全用它的抽样值来代表,而不会丢掉如何信息。这就是抽样定理的意义。

    总结:如果序列是通过对模拟信号采样得到的,有关系:x(n)=xa(nT),即序列值对于对模拟信号的采样值,或者说对于采样信号在t=nT时的幅度。

    例:2724041_220.gif= sin(Wt),理想抽样后,

    x(n)= 2724041_224.gif=sin(WnT)= sin(nω0)

    ∴ω0=WT数字域频率与模拟角频率之间的关系。

    ω0=WT=W/fs=2pf/fs  2. 4时域离散信号的傅立叶变换与模拟信号傅立叶变换之间的关系

    连续信号2724041_220.gif的傅立叶变换及反变换公式如下:

    2724041_225.gif

    2724041_220.gif=2724041_226.gif

    理想抽样后的抽样信号为2724041_190.gif

    2724041_190.gif= xa(t) p(t) =2724041_208.gif

    则抽样信号的傅立叶变换

    2724041_213.gif = 2724041_205.gif

    离散时间信号x(n)=xa(nT),x(n)的傅立叶变换为:2724041_1.gif(2.2.1)

    抽样信号的傅立叶变换2724041_213.gif2724041_227.gif有什么关系?

    可以证明,2724041_227.gif也可写成:

    2724041_227.gif=2724041_228.gif,对照2724041_213.gif

    2724041_213.gif = 2724041_205.gif,都是2724041_229.gif以周期2724041_197.gif进行周期延拓。

    2724041_227.gif时,以w为横轴,以周期2724041_230.gif进行周期延拓。

    2724041_213.gif时,以W为横轴,以周期2724041_197.gif进行周期延拓。

    坐标轴之间的对应关系如下图所示。

    2724041_231.jpg

    在一些文献中经常使用归一化频率2724041_232.gif,因为2724041_233.gif2724041_234.gif2724041_235.gif都是无量纲量,刻度是一样的。

    所以数字频率0、2p处是低频,p附近代表高频。

    当抽样频率是信号最高频率2724041_236.gif4倍时,最高频率2724041_236.gif所对应的数字频率为2724041_237.gif

    展开全文
  • 西华师范大学 China West Normal University 1 第七章 信号频域分析及 MATLAB 实现 7.1 周期信号的傅利叶级数与信号的频谱 7.2 周期信号的频谱分析及 MATLAB 实现 7.3 用 MATLAB 分析典型周期信号的频谱 西华师范...
  • 利用Matlab信号进行频域分析 6.1离散时间周期信号离散时间傅里叶级数 基本周期为N基频为0=2/N周期信号x[n]的离散时间傅里叶级数为 其中 在matlab中可以用fft和ifft来求解离散时间傅里叶级数长度为N矢量x可以...
  • 尚戈继:连续时间滤波器Matlab仿真实例​zhuanlan.zhihu....给出了代码本篇目录:待处理信号离散时间滤波器滤波器指标IIR滤波器设计FIR滤波器设计FIR和IIR效果对比IIR和FIR对比完整代码参考待处理信号与连续时间...

    abccb4b2bf03d94ef9c0d6a5242cd759.png
    尚戈继:连续时间滤波器Matlab仿真实例zhuanlan.zhihu.com
    0845705a02bc87b9f174b91ffe1bc936.png
    本文是这篇文章的后续,讨论了两类离散时间滤波器(无限冲激响应IIR和通过Kaiser窗函数法得到的有限冲激响应FIR)在Matlab中的实现。进行了对比分析。给出了代码

    本篇目录:

    • 待处理信号
    • 离散时间滤波器
      • 滤波器指标
      • IIR滤波器设计
      • FIR滤波器设计
      • FIR和IIR的效果对比
    • IIR和FIR的对比
    • 完整代码
    • 参考

    待处理信号

    与连续时间情况一样,待处理的信号是gamme.m函数生成的一段依次从Do、 Re、Mi、...一直到升Do的音频,每个频率持续一秒钟。

    function [gam,t] = gamme(duree,fe)
    
    % [do re mi fa sol la si doo]
    freqnotes = [262 294 330 349 392 440 494 523];
    t = 0:1/fe:duree;
    
    gam = [];
    for note=1:8, 
        gam = [gam, sin(2*pi*freqnotes(note)*t)];
    end
    N = length(gam);
    t = (0:N-1)/fe;

    调用函数,生成音频,并使用matlab的soudnsc函数试听

    clear all
    fe = 8192 ; %采样频率
    duree = 1; %每个音符持续1秒
    [sig,t] = gamme(duree,fe);
    soundsc(sig, fe); % listen

    对原始信号进行频率分析

    fmin = 250;
    fmax = 550;
    
    f = [0:length(sig)-1]*(fe/length(sig));
    S = fft(sig)/fe;
    
    figure()
    plot(f, abs(S));
    axis([fmin fmax 0 0.5])

    e7e0e3d12a90f6920904edd6f884149b.png
    八个音,每个的幅值都是一样的,是因为它们在信号中的比例(时长)相同

    离散时间滤波器

    滤波器指标

    与连续时间情况不同,此时我们希望除去FA这个音,此时的滤波器指标

    7bafdc33b080cf079f19db08f8952643.png

    类型:陷波滤波器

    通带容许的最小增益:-1dB

    阻带增益:-40dB

    过渡带中心频率:340Hz和360Hz(【Mi和Fa的中点】和【Fa和Sol的中点】)

    过渡带宽度:10Hz

    d6e3a25fe2d0b13232bca188fe6d212c.png
    陷波滤波器参数示意图,注意不同滤波器的通带和阻带是可能存在抖动的,只是图中没有体现出来

    IIR滤波器设计

    0、IIR滤波器介绍

    这里引用一段话[1]

    历史上,随着数字信号处理领域的兴起,离散时间IIR滤波器的设计依赖于将连续时间滤波器变换成满足预定指标的离散时间滤波器的方法。该方法过去是,并且现在仍然是一种合理的方法, 其理由如下:
    •连续时间IIR滤波器的设计技巧十分成熟,并已取得许多有用的成果.因此可以方便地利用这些为连续时间滤波器推导出的设计方法。
    •许多有用的连续时间滤波器设计方法有比较简单的完整设计公式。因此,以这种标准的连续时间IIR滤波器设计公式为基础的离散时间IIR滤波器的设计方法实现起来十分简单。
    •当把完全适用于连续时间IIR滤波器的标准逼近方法直接用于离散时间IIR滤波器时,并不能得出简单的完整设计公式,因为离散时间滤波器的频率响应是周期的,而连续时间滤波器则不然。

    1、常量定义

    fcb = 340/(fe/2); %低频段的
    fch = 360/(fe/2); %高频段的
    deltaf = 10/(fe/2);
    Wp = [fcb-deltaf/2,fch+deltaf/2]; %pass 通带的两个边界
    Ws = [fcb+deltaf/2,fch-deltaf/2]; %stop 阻带的两个边界
    
    Ap = 1;
    Aa = 40;
    • 这里Wp和Ws的定义见Cheby的文档 Chebyshev Type II filter order

    3a7e9e95340312c3c6c1434cd5006d4d.png
    • 在离散的情况下频率需要做出对应的变换,即 离散频率=实际频率/(fe/2),其中fe代表信号采样频率

    2、计算滤波器阶数

    • 这里我们使用的连续滤波器原型是Chebyshev Type II filter,在通带无抖动,在阻带有抖动
    • 当然也可以选用其他类型的滤波器,方法类似,在连续时间滤波器Matlab仿真实例中介绍过了
    [nd wnd]= cheb2ord(Wp,Ws,Ap,Aa) %这里不写‘s’就是离散的 %命名中的‘d’代表discrete
    • 输出

    5026ef1bdf8ffbf8c7717a2f8af0245f.png
    只需要4阶就能实现,还是比较容易的
    • 注意这里得到的两个wnd是通带和过度带的交界位置频率。

    3、计算滤波器参数

    [numd, dend] = cheby2(nd,Aa,wnd,'stop')

    可以画出这个滤波器的频域响应图

    f = linspace(200,550,1000);
    freqz(numd,dend,f,fe)

    d819647c2ab87f4a68338b75a0072d1a.png
    • 得到的滤波器很好地满足了设计指标,在阻带的增益小于-40(标的那条线),在增带的增益大于-1dB(图中看不太出来)
    • 也可以看到IIR的缺点:相位是非线性的
    • 滤波效果放在下下节,和FIR进行对比

    FIR滤波器设计

    0、FIR滤波器介绍

    通常设计IIR滤波器采用的方法都是由连续时间系统到离散时间IIR系统的变换演变而来的。与之对应,离散FIR滤波器的设计方法以直接逼近所需离散时间系统的頻率响应或脉冲响应为基础。

    设计FIR滤波器的最简单的方法称为窗函数法,核心思想在于:根据预期需求,在频域内设计一个“窗“函数--->得到该”窗“的脉冲响应,(非因果和无限长的)--->截断该理想脉冲响应,得到预期的FIR滤波器。

    cf15bb9614954dbbcf30fa2f39e6a3b6.png

    最容易想到的窗函数自然就是分段恒定的窗,如上图所示,但是,这种窗会引起吉布斯现象,即截断得到的FIR滤波器会随着有限响应长度(即截断的长度)增加,非一致地收敛到窗函数,而这种非一致收敛并不是我们想要的,具体来讲,这时滤波器并不能有效地滤掉不想要的频率成分,无论增加多少有限响应的长度。

    因此,实际中并不会使用上图所示的窗,而是别的窗,如Bartlett、Hanning、Hamming、Blackman等,这些窗口的两端平滑地缩小到零,吉布斯现象会被减弱,但代价是主瓣加宽,过渡带的带宽增大,因此为了达到相同的性能指标,需要更高的阶数(即长度)。

    3ec0f650dcb6ce446968b8db9d290ffe.png
    常用的窗函数

    这里使用的窗函数是Kaiser窗,与上图所示的窗不同,Kaiser窗有两个自由度,它是一种近似最佳的窗函数。

    1、常量定义

    这里不用做向离散的频率变换

    fcb = 340;
    fch = 360;
    deltaf = 10;
    Wp = [fcb-deltaf/2,fch+deltaf/2];
    Ws = [fcb+deltaf/2,fch-deltaf/2];
    
    Ap = 1;
    Aa = 40;

    2、滤波器系数计算

    [n,wc,beta] = kaiserord([fcb-deltaf/2 fcb+deltaf/2 fch-deltaf/2 fch+deltaf/2]...
        ,[1 0 1],[1-10^(-Ap/20) 10^(-(Aa)/20) 1-10^(-Ap/20)],fe)

    [n,Wn,beta,ftype] = kaiserord(f,a,dev,fs)

    f:频率区间,a:不同区间的预期增益,[1 0 1]代表着[通 阻 通]

    dev:不同区间内相较预期增益的容许偏差 我们的性能指标是以分贝的形式给的,这里要转换成小数

    fe:采样频率

    5224768d45786e09f9652f85ac90dbee.png
    计算结果

    3、滤波器

    使用firl截取窗函数的脉冲响应,获得FIR滤波器

    h = fir1(n,wc,'stop',kaiser(n+1,beta));
    • Kaiser窗函数的设置保证了得到的脉冲是因果并有广义线性相位特性的。因此滤波器的传递函数的分子系数num=h,分母系数den=1

    同样的,画出该滤波器的频率响应

    f = linspace(200,550,1000);
    freqz(h,1,f,fe) %num = h den = 1

    526260fb5975d3787d2495600869789d.png
    Kaiser窗函数+截取法得到的FIR的频域相应

    FIR和IIR的效果对比

    频域响应上

    figure()
    hold on 
    FIR = freqz(h,1,f,fe); %FIR
    IIR = freqz(numd,dend,f,fe); %IIR
    
    hold on
    grid on 
    plot(f,20*log10(abs(IIR)))
    plot(f,20*log10(abs(FIR)))
    legend('IIR','FIR')

    f4e4357b6c591988cae81af84a425007.png
    可以看到IIR在阻带严格的小于-40dB,而FIR在阻带有一定的误差(红色偏高一点)

    Tips: FIR的误差大小约在3dB左右,因此,如果想让阻带的增益更加精确地小于-40dB,一个常用的技巧是在计算Kaiser窗函数系数时预先补偿这3dB,即

    [n,wc,beta] = kaiserord([fcb-deltaf/2 fcb+deltaf/2 fch-deltaf/2 fch+deltaf/2]...
        ,[1 0 1],[1-10^(-Ap/20) 10^(-(Aa+3)/20) 1-10^(-Ap/20)],fe)% 原本是10^(-(Aa)/20),改成10^(-(Aa+3)/20)

    就能得到

    1a61b96f4a05f7212c203111140855cb.png

    滤波效果上

    filter:调用离散时间滤波器处理信号

    sigf = filter(numd,dend,sig); %IIR
    sigf2 = filter(h,1,sig); %FIR %num = h den = 1

    然后看输出的信号的频域成分

    • FFTplot函数的定义详见连续时间滤波器Matlab仿真实例文末,其实很简单,就是画出一个信号的傅里叶变换图,并放大到特定区间。
    figure()
    subplot(3,1,1)
    FFTplot(sig, fe,200, 550) 
    title('raw')
    
    subplot(3,1,2)
    FFTplot(sigf, fe,200, 550)
    title('IIR')
    
    subplot(3,1,3)
    FFTplot(sigf2, fe,200, 550)
    title('FIR')
     

    4053c460628c8a9e52c14dfc9bce92c1.png

    差不多,但是输出信号在时域上则有不同

    figure()
    subplot(3,1,1)
    plot(t, sig)
    title('raw')
    ylim([-2 2]);
    
    subplot(3,1,2)
    plot(t, sigf)
    title('IIR')
    ylim([-2 2]);
    
    subplot(3,1,3)
    plot(t, sigf2)
    title('FIR')
    ylim([-2 2]);

    e9b4404e464fa5c5491e9e045871f5a3.png
    • 可以看到FIR有时延的效果。这源于它的广义非线性相位特性
    • IIR并没有时延的效果,IIR的相位时非线性的,得到的信号相位信息会被干扰,因此IIR适用于那些对相位信息不敏感的信号上。
    • IIR对信号相位信息的干扰,体现在
      • IIR处理得到的信号的”毛刺“更多一点(第4、5、6秒处)
      • IIR的通带和阻带的边界的滤波效果更差一点(第3秒处)
    • 这种对相位信息的干扰,对于我们这里处理的相对比较简单的信号影响十分小。
    • 可以试听一下滤波得到的声音
    %不要三行命令同时运行,一次只运行一个
    %如果电脑有杜比全景声(Dolby Atmos)的话,暂时关掉再播,否则听到的声音有失真
    soundsc(sig, fe); % RAW listen
    % soundsc(sigf, fe); % IIR listen
    % soundsc(sigf2, fe); % FIR listen

    IIR和FIR的对比

    标准FIRIIR
    相位可控不可控
    复杂度相对较高
    稳定性恒稳定不恒稳定
    阶数或参数个数
    精准程度普通
    • 值得注意的是,虽然FIR在直观上是比IIR更加简单实现的,但是它的精确程度比IIR差,因此,对于同样的滤波器设计指标,FIR滤波器所要求的阶数可能比IIR滤波器高5-10倍,结果,成本较高,信号延时也较大;如果按线性相位要求来说,则IIR滤波器就必须加全通网络进行相位校正,同样要大大增加滤波器的阶数和复杂性。而FIR滤波器却可以得到严格的线性相位。[2]
    • IIR滤波器必须采用递归结构,极点位置必须在单位圆内,否则系统将不稳定。另外,在这种结构中,由于运算过程中对序列的舍入处理,这种有限字长效应有时会引入寄生振荡。相反,FIR滤波器主要采用非递归结构,不论在理论上还是在实际的有限精度运算中都不存在稳定性问题,运算误差也较小。此外,FIR滤波器可以采用快速付里叶变换算法,在相同阶数的条件下,运算速度可以快得多。[2]

    以下三段来自[1]

    • 选择IIR或是FIR滤波器取决于每种类型滤波器的优点在设计问题中的重要性。例如,IIR 滤波器具有可以用完整的设计公式来设计各种选频滤波器的优点。这就是说,一旦选定了选用哪种已知的逼近方法(即巴特沃思、切比雪夫或椭圆逼近),则可以直接把技术指标代入一组设计方程来计算满足技术条件的滤波器的阶次,并得出数字滤波器的系数(或极点和零点)。这种IIR滤波器的非迭代计算程序。这些方法只限于设计选频滤波器,并只允许用于规定了幅度响应的场合。如果要得到其他形状的幅度响应,或需要逼近预定的相位响应或群延迟响应。则需要用算法设计法。
    • 与此相反,FIR滤波器可以有精确的(广义)线性相位。但是对于FIR滤波器不存在完整的设计方程。虽然可以直接用窗函数数法,但为了满足预定的技术指标,有可能需要做一些迭代。与窗困数法相比。Parks-Mcclellan算法可以得出较低阶的滤波器,并且这两种方法的墟波器设计程序都很容易得到。而且。窗函数法和大多数算法设计法都有可能逼近较为任意的频率响应特性,但所遇到的困难要比在低通浦波器设计中遇到的稍大一些。。此外,FIR滤波器的设计问题要比IIR设计问题有更多的可控之处,因为对于FIR滤波器有适用于各种实际情况的最佳理论。
    • 最后,在实现数字滤波器时还要考虑成本问题。通常将硬件的复杂性。芯片的面积或计算速度等作为衡量成本问题的因素。这些因素或多或少地直接与满足给定指标所需的滤波器阶次有关。在实际应用中,多相位实现的功效并没有表现出来,一般说来用IIR滤波器就能最有效地满足给定的幅度响应技术指标。但是在许多情况下,FIR滤波器的线性相位与它所带来的额外成本相比是非常值得的。

    简单地说,FIR设计比较简单,IIR更加灵活,但也更加复杂。

    完整代码

    %信号生成
    clear all
    fe = 8192 ;
    duree = 1;
    [sig,t] = gamme(duree,fe);
    soundsc(sig, fe); % listen
    fmin = 250;
    fmax = 550;
    
    %常量定义
    fcb = 340/(fe/2);
    fch = 360/(fe/2);
    deltaf = 10/(fe/2);
    Wp = [fcb-deltaf/2,fch+deltaf/2];
    Ws = [fcb+deltaf/2,fch-deltaf/2];
    
    Ap = 1;
    Aa = 40;
    
    
    %%
    %IIR 
    %滤波器参数
    
    [nd wnd]= cheb2ord(Wp,Ws,Ap,Aa)
    [numd, dend] = cheby2(nd,Aa,wnd,'stop')
    
    f = linspace(200,550,1000);
    figure()
    freqz(numd,dend,f,fe)
    grid on
    %滤波
    sigf = filter(numd,dend,sig);
    
    %%
    %FIR 
    %常量定义
    fcb = 340;
    fch = 360;
    deltaf = 10;
    Wp = [fcb-deltaf/2,fch+deltaf/2];
    Ws = [fcb+deltaf/2,fch-deltaf/2];
    
    Ap = 1;
    Aa = 40;
    
    %滤波器参数
    E_comp = 3; % RIF最后在阻带会有3dB左右的误差,可以通过预先补偿的方式消除
    [n,wc,beta] = kaiserord([fcb-deltaf/2 fcb+deltaf/2 fch-deltaf/2 fch+deltaf/2]...
        ,[1 0 1],[1-10^(-Ap/20) 10^(-(Aa+3)/20) 1-10^(-Ap/20)],fe)
    % dev is a vector the same size as a that specifies the maximum allowable error...
    %     or deviation between the frequency response of the output filter and its...
    %     desired amplitude, for each band. The entries in dev specify the passband...
    %      ripple and the stopband attenuation.
    h = fir1(n,wc,'stop',kaiser(n+1,beta));
    freqz(h,1,f,fe);
    
    sigf2 = filter(h,1,sig);
    
    %%
    %画图
    
    figure()
    hold on
    FIR = freqz(h,1,f,fe);
    IIR = freqz(numd,dend,f,fe);
    
    hold on
    grid on
    plot(f,20*log10(abs(IIR)))
    plot(f,20*log10(abs(FIR)))
    legend('IIR','FIR')
    
    figure()
    subplot(3,1,1)
    FFTplot(sig, fe,200, 550)
    title('raw')
    
    subplot(3,1,2)
    FFTplot(sigf, fe,200, 550)
    title('IIR')
    
    subplot(3,1,3)
    FFTplot(sigf2, fe,200, 550)
    title('FIR')
    
    figure()
    subplot(3,1,1)
    plot(t, sig)
    title('raw')
    ylim([-2 2]);
    
    subplot(3,1,2)
    plot(t, sigf)
    title('IIR')
    ylim([-2 2]);
    
    subplot(3,1,3)
    plot(t, sigf2)
    title('FIR')
    ylim([-2 2]);
    
    %%
    %试听
    % soundsc(sig, fe); % RAW listen
    % soundsc(sigf, fe); % IIR listen
    % soundsc(sigf2, fe); % FIR listen

    参考

    1. ^ab离散时间数字信号处理,奥本海姆 https://book.douban.com/subject/26307919/
    2. ^abhttps://rf.eefocus.com/article/id-332067
    展开全文
  • spm_id_from=pageDriver 第二章 信号函数与标准信号 ...连续正弦波要进行离散化,才能进而进行数字信号处理,将连续的时间t用n倍的时间间隔表示 改错:dt=1/11025 matlab中生成三角波和锯齿

    【学习笔记】matlab进行数字信号处理(一)生成信号及信号的时域频域分析
    【学习笔记】matlab进行数字信号处理(二)信号的相关分析及幅值分析
    【学习笔记】matlab进行数字信号处理(三)数字滤波技术
    【学习笔记】matlab进行数字信号处理(四)信号的时频域分析

    b站视频地址:https://www.bilibili.com/video/BV18E411f7ZQ?p=16&spm_id_from=pageDriver

    第二章 信号函数与标准信号

    2.1 概述

    信号发生器产生信号,可用信号函数生成正弦波、方波、三角波、锯齿波、白噪声、脉冲信号、阶跃信号、斜波信号、加速度信号

    2.2信号函数和标准信号

    在这里插入图片描述
    连续的正弦波要进行离散化,才能进而的进行数字信号处理,将连续的时间t用n倍的时间间隔表示
    在这里插入图片描述
    改错:dt=1/11025
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    matlab中生成三角波和锯齿波用同一个函数sawtooth
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    2.3 信号发生器设计

    在这里插入图片描述
    在这里插入图片描述
    b站视频P20:讲解如何生成信号发生器和电子琴(包含GUI设计)

    第三章 信号的时域分析

    3.1概述

    在这里插入图片描述

    3.2信号波形参数识别

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    画一条零线,可求出信号的周期和初相位
    相邻两过零点时间差就是周期,第一个过零点位置与周期的比例可算出初相位
    过零检测法求周期和相位很常用
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    3.3 信号的数字微分/积分

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述在这里插入图片描述

    第四章 信号的频谱分析

    在这里插入图片描述

    4.1频谱分析的概念

    在这里插入图片描述
    在这里插入图片描述
    可以从频谱图看出一个信号的频率构成,和不同频率不分的强弱(幅值)
    从波形,可以看出信号随时间变化的幅值强弱
    与波形相比,频谱可以更直观的
    波形:幅值随时间的变化情况;频谱:幅值随频率的变化情况
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    与波形相比,频谱的抗干扰能力更强
    在这里插入图片描述
    频谱分析分解出的正弦波,往往能找到明确的物理意义,比如齿轮转动的频谱,每一个谱线都能对应到一个机械零件上去

    4.2 周期信号的频谱分析(FFT代码和作业)

    在这里插入图片描述
    在这里插入图片描述
    任意两个波形相乘,积分等于0说明两个函数正交
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    功率谱是幅频谱的平方,反映了每个频率分量的能量大小
    在这里插入图片描述
    周期信号的频谱图:
    谐波性是指频率成分一定是基频的倍数,如基频是f0,则谱线一定智能出现在f0的整数倍上
    收敛性是指谐波次数越高,幅值越小
    在这里插入图片描述
    在这里插入图片描述
    fft使用条件:数据长度必须是2的次方
    在这里插入图片描述
    合成信号——FFT——实频谱和虚频谱——幅值谱和相位谱
    在这里插入图片描述
    【视频P27,11分位置】注意:使用FFT之后,并不需要显示负频谱的部分,因为正频谱部分就可以显示出有用的信息,负频谱部分的信息是冗余的
    在这里插入图片描述
    有时功率谱上某些频谱分量的幅值显示的不明显,小能量信号会被大能量信号掩盖,所以通常采用对数功率谱,能将小的信号放大,大的信号压缩,能在一张图上表现出不同大小的频率成分
    在这里插入图片描述
    在这里插入图片描述

    4.3 数字信号的频谱计算方法

    在这里插入图片描述

    在这里插入图片描述
    截取后再延拓的信号会在连接点处产生跳变,在波形上显示的是产生跳变,在频谱上产生的是出现能量泄露
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    4.4 FFT中的能量泄露和栅栏效应

    在这里插入图片描述
    快速傅立叶变换和离散傅立叶变换的计算结果是相同的,但是速度快很多
    在这里插入图片描述
    在这里插入图片描述
    栅栏效应生成的原因:为了提高计算效率,在使用FFT计算频谱时,在0到二分之一的采样频率区间内,不是每个频率点都计算,而是按照一定的频率间隔抽样计算,频率间隔等于采样频率除以数据长度
    在取样的位置和谱峰的位置不重叠时,便会带来一个谱峰误差,这种误差引起的反应,叫做栅栏效应。包含两个误差,一个是幅值高度的误差,一个是最大频率的误差。
    在这里插入图片描述
    谱峰越尖锐,产生误差的可能性就越大
    在这里插入图片描述
    在这里插入图片描述
    由于截断带来的能量泄露,这时会减小栅栏效应的误差,可以通过控制截断函数的形状,来调整能量泄露形状,
    左侧图为矩形窗能量泄露函数,以中心频率为主的成份叫做主瓣,比较窄比较尖,越尖锐造成的栅栏效应的误差就越大,还有其他很高的旁瓣,但理想情况下只希望读出主瓣处的能量泄露,但这种理想的窗函数并不存在,实际上只能去逼近它。
    在这里插入图片描述

    直接截断,然后周期延拓,如果截断的不是整周期的位置,会产生一个跳变,此时加一个汉明窗,用这样的窗函数和原来的信号相乘,这样会将原信号在接头处的幅值被压缩成0,以这种方式进行周期延拓,无论原来的信号怎样,在接头处都不存在跳变(从波形上看)
    (从频谱上看)相当于增宽主瓣,压缩旁瓣,因此可以通过这种方式抑制旁瓣的能量泄露,同时提高主瓣的宽度,使栅栏效应的误差减小
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    汉宁窗和矩形窗相比,主瓣变的更宽,而旁瓣变得更低
    在这里插入图片描述
    这是汉宁窗对波形进行截断的情况,截断后要进行修正,窗函数和原信号相乘再进行截断延拓,幅值修正后幅值增大了一倍,这样幅值增大后可以弥补两端衰减造成的能量损失
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    和矩形窗相比,主瓣宽一些,旁瓣窄一些
    在这里插入图片描述
    主瓣更宽一些,旁瓣更窄一些
    在这里插入图片描述
    在这里插入图片描述
    平顶窗,主瓣相当宽
    在这里插入图片描述
    在这里插入图片描述
    能量修正z=2*w1.*x,然后再进行FFT
    加窗能克服能量泄露和栅栏效应带来的误差
    在这里插入图片描述

    4.5 非周期信号的频谱分析

    在这里插入图片描述
    工程上大量的信号都是非周期信号,在数字领域,由于截断,不管原信号是不是周期信号,都会转换成周期信号
    在这里插入图片描述
    1.如果数据长度不够长,可以通过补0的方式提高频谱的频率分辨率
    eg:假如原来数据长度是1024,再补1024个0,频谱分析时,精度会高一倍
    2.细化:如果想看清某一部分的频谱,使用ZOOM-FFT可以将观测的视角集中在一个频率
    3.主要针对受噪声干扰频谱,计算时谱当中有很多的随机干扰,此时将频谱多次累加再除以N,这样干扰也会减小到N分之一
    4.用于修正栅栏效应误差
    5.传感器输出的都是实信号,但是matlab里的FFT输入的是复信号,因此针对虚部的计算浪费了,采用实信号FFT计算技术,可以将计算速度提高一倍

    4.6频谱分析应用

    在这里插入图片描述
    在这里插入图片描述
    结果会发现,信号波形的幅值是1,倒是频谱的幅值不是1,产生差异的原因就是能量泄露和栅栏效应

    大作业:声音信号采集和频谱分析程序设计

    首先通过计算机上的脉冲,采集一段声音信号,然后进行FFT变换,再画出信号频谱
    1.新建GUI文件
    2.安放控件,修改控件属性
    3.使用定时器,实现连续采样
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    展开全文
  • 这篇博文中使用的模拟信号为上篇博文:【 MATLAB 】使用 MATLAB 实现模拟信号的近似及其连续傅里叶变换 中使用的模拟信号: 为了研究在频域数量上的采样效果,对该信号使用两种不同的采样频率采样。 a. 在 fs =...
  • 数字信号处理DSP基本原理之一:采样时间信号的频谱具有周期性,且周期与采样率相等。这个结论可以帮助我们在大脑中构建这样的一幅图像:时域的动态采样,相当于在频域进行周期延拓,延拓的周期与时域采样率相等。教...
  • 信号与系统离散频域分析(DFT) 一、实验目的: 1、掌握离散时间系统DFT的MATLAB实现; 2、熟悉DTFT和DFT之间关系。 3、了解信号不同变形DFT与原信号DFT之间关系
  • 信号重建是将离散信号x[k]转换为连续时间信号x(t)。 非周期离散信号的频谱是连续的周期谱。计算机在分析离散信号的频谱时,必须将其连续频谱离散化。频域抽样定理给出了连续频谱抽样过程中信号不失真的约束条件。 ...
  • 奥本海姆所著《信号与系统》(刘树棠译版)中关于Parseval定理的描述如下:信号的能量既可以按每单位时间的能量在整个时间内积分出来,也可以按每单位频率的能量在整个频率范围内积分而得到。简单地讲,信号从时域...
  • (本人为本文原作者,转载请标明出处)上文请转:DBinary:傅里叶变换推导详解​zhuanlan.zhihu.com其中,式3.51及3.52被称为傅里叶变换对,其中3.52为正变换,表示离散时间信号变换为频域复信号,式3.51位逆变换,表示由频域...
  • DFT源头,是连续傅里叶变换,用于将连续时间信号x(t)转换成连续频域信号X(f)。但是,连续傅里叶变换不适合计算机上应用,所以工程师们就发明了离散傅里叶变换(DFT)。其中,x(n)是离散输入序列,X...
  • 3.深刻理解离散时间系统系统函数在分析离散系统时域特性和频率特性中重要作用及意义。 二、实验内容 1、已知系统y(k)-y(k-1)-2y(k-2)=f(k),输入f(k)=(-1)^k ε(k),起始状态为y(-1)=0, y(-2)=1/6。 (1)...
  • 连续时间信号的频域分析及matlab实现 离散时间系统的频域分析及matlab实现 傅里叶变换 零极点分析 全部附有编译通过的代码
  • 之前文章信号频域分析方法理解(频谱、能量谱、功率谱、倒频谱、小波分析)中提到了离散小波分解例子,其参考代码如下:t_s = 0.005; %采样周期 t_start = 0.001; %起始时间 t_end = 10; %结束时间 t = t_...
  • 离散时间信号——序列——可以用图形来表示。 按信号特点不同,信号可表示成一个或几个独立变量函数。例如,图像信号就是空间位置(二元变量)亮度函数。一维变量可以是时间,也可以是其他参量,习惯上将其...
  • 掌握离散时间信号与系统时域分析方法。2.掌握序列傅氏变换计算机实现方法,利用序列傅氏变换对离散信号、系统及系统响应进行频域分析。3.熟悉理想采样性质,了解信号采样前后频谱变化,加深对采样定理...
  • 目 录 绪论 1 离散时间信号和系统分析 1.1 离散时间信号产生与运算 2 1.2 离散时间系统的时域分析 9 1.3 离散时间系统的频域分析 13 1.4 离散时间系统频响的零极点确定 14 快速傅立叶变换的应用 2.1 FFT的计算 17 ...
  • DFT和FFT的区别1 连续时间信号频域分析2 通过离散时间信号的Z变换表达式X(z)直接写出时域离散信号(序列)x(n)的方法3 部分分式法的MATLAB实现(求X(z)的部分展开式)4 稳定系统5 求频响特性(系统函数H与对应的频点w)6 ...
  • 当我们处理一段信号时,如图1所示,往往需要通过消除噪声,放大有效信号的部分的处理,此时可以通过时域和频域的方式来滤除噪声。FIR与IIR的区别在于,FIR具有严格的线性相位特性,但运算时间过长。而IIR滤波器则...
  • matlab在DSP中的应用(四)---时域抽样与信号的重建

    万次阅读 多人点赞 2017-04-05 12:33:43
    一、实验目的(1)掌握用MATLAB语言进行离散时间傅里叶变换和逆变换方法。...二、实验原理1.DTFT离散时间傅里叶变换(DTFT)是指信号在时域上为离散,而在频域上则是连续。如果离散时间非周期信号为x(n),
  • 之前文章信号频域分析方法理解(频谱、能量谱、功率谱、倒频谱、小波分析)中提到了离散小波分解例子,其参考代码如下: t_s = 0.005; %采样周期 t_start = 0.001; %起始时间 t_end = 10; %结束时间 t = t_...
  • 摘? 要 介绍了利用TLAB信号处理工具箱进行FI滤波器设计的三种方法... 数字滤波器是一种用来过滤时间离散信号的数字系统,通过对抽样数据进行数学处理来达到频域滤波的目的根据其单位冲激响应函数的时域特性可分为两类:
  • 傅里叶变换可以将时域信号转变成频域,通过分析频谱了解信号的组成。网上有大量介绍傅里叶变换的好文章,感兴趣的小伙伴可以自行查阅!什么是时域和频域呢?简单的理解是:时域的横轴为时间,反映信号随时间的变化,...
  • 连续时间 LTI 系统的时域分析 实验四 傅里叶变换系统的频域分析 实验五 信号抽样与恢复 实验六 信号与系统复频域分析 实验一 基本信号在 MATLAB 中的表示和运算 一实验目的 1 学会用 MA TLAB 表示常用连续信号的方法
  • 摘 要 介绍了利用MATLAB信号处理工具箱进行FIR滤波器设计的三种方法:程序设计法、FDATool设计法和SPTool设计法,... 数字滤波器是一种用来过滤时间离散信号的数字系统,通过对抽样数据进行数学处理来达到频域滤波的
  • matlab开发-频率域零填充重新采样Interpolation。离散时间信号的频域(基于FFT)重采样

空空如也

空空如也

1 2 3 4 5
收藏数 86
精华内容 34
关键字:

matlab离散时间信号的频域

matlab 订阅