matlab编写信号的变换

前言

信号与系统的分析方法中,除了时域分析方法外,还有变换域分析的方法。连续时间信号与系统的变换域分析方法主要是傅里叶变换和拉普拉斯变换。离散时间信号的Z变换是分析线性时不变离散时间系统问题的重要工具,在数字信号处理、计算机控制系统等领域有着广泛的应用。
学习目标:

• (1)了解、熟悉Z变换的概念与性质;
• (2)理解Z反变换的相关内容;
• (3)掌握离散系统中的Z域描述方法;
• (4)了解、熟悉傅里叶级数与变换;
• (5)理解离散傅里叶变换及其性质;
• (6)实现频率域采样和快速傅里叶变换;
• (7)熟悉实现离散余弦变换.ChirpZ变换和Gabor函数。

1.Z变换

这里就不在赘述Z变换的性质和定理。

MATLAB的符号数学工具箱提供了计算Z变换的函数ztrans()和Z反变换的函数iztrans(),其调用形式为
F= ztrans(f),f= iztrans(F)
其中,右端的f和F分别为时域表示式和Z域表示式的符号表示,可应用函数sym来实现,其调用格式为S=sym(A)，
在MATLAB中,留数法求Z反变换可以使用函数residuez实现，调用格式为[R P K] = residuez(B,A);
其中,B和A分别为X(z)的多项式中分子多项式和分母多项式的系数向量;返回值R为留数向量,P为极点向量,二者均为列向量;返回值K为直接项系数,仅在分子多项式最高次幂大于等于分母多项式最高次幂时存在,否则,返回值为空。

例1 求f(n)=sin(ak)u(k)的Z变换和F(z)=z/(z-3)^2的Z反变换。

syms f(k) a k
f= sin(a*k);
F=ztrans(f)
syms F(z) 
F= z/(z-3)^2;
f= iztrans(F)

结果如下:

例2 求X(z)=(1+0.4z^-1)^2/ (1+0.8z^-1)^2* (1-0.5z^-1)^2*(1+0.1z ^-1), |z|>0.8 的Z反变换。

clear all;close all;clc;
B=poly([-0.4 -0.4]);
A=poly([ -0.8 -0.8 0.5 0.5 -0.1]);
[R P K]= residuez(B, A);
R= R'
P= P'
K

结果如下：

2.离散系统函数的频域分析

在MATLAB中，函数freqz用于求离散时间系统频响特性，该函数的调用方法如下:
[H,w] = freqz(B, A,N)
[H,w] = freqz(B, A, N, ‘whole’)
其中,B与A分别表示H(z)的分子与分母多项式的系数向量; N为正整数,默认值为512;返回值w包含[0,pi]范围内的N个频率等分点;返回值H则是离散时间系统频率响应H(e^jw)在0~π范围内N的频率处对应的值。

例3 绘制如下系统的频响曲线:H(z) =z- 1.3/z

B=[1 -1.3];
A =[1 0];
[H,w]= freqz(B, A, 400, 'whole');
Hf = abs(H);
Hx= angle(H);
clf
subplot(121)
plot(w, Hf)
title('离散系统幅频特性曲线')
xlabel( '频率');ylabel('幅度')
grid on
subplot(122 )
plot(w, Hx)
xlabel( '频率');ylabel('幅度')
grid on
title('离散系统相频特性曲线')

例4 绘制离散系统的幅频响应和相频响应

clear all;close all;clc;
w=(-4*pi:0.001:4*pi) + eps;
X=1./(1-0.6* exp( -j*w));
subplot(211),
plot(w/pi, abs(X),'LineWidth',2);
xlabel( '\omega/\pi');
ylabel('|H(e^j^\omega)|');
title('幅频响应');
axis([-3.2 3.2 0.5 2.2]);
grid;
subplot(212),
plot(w/pi, angle(X), 'LineWidth' ,2);
xlabel('\omega/\pi');
ylabel('\theta(\omega)');
title( '相频响应');
axis([ -3.2 3.2 -0.6 0.6]);grid;
set(gcf, 'color', 'w');

3.离散系统函数零点分析

在MATLAB中，系统函数的零极点就可以通过函数roots得到，也可以借助函数tf2zp得到,tf2zp的语句格式为
[Z, P , K] = tf2zp(B, A)
其中,B与A分别表示为H(z)的分子与分母多项式的系数向量。上式的作用是将H(z)
的有理分式表示为转换为零极点增益形式,即
H(z) = k (z-z1)(z-z2)…(z-zm ) /(z- p1)(z-p2)…(z-pn)
zplane函数用于绘制H(z)的零极点图,该函数的调用格式为
zplane(z,p)
绘制出列向量z中的零点(以符号“o”表示)和列向量p中的极点(以符号“X"表示),以及参考单位圆,在多阶零点和极点的右上角标出其阶数。如果z和p为矩阵,则会以不同颜色绘出z和p各列中的零点和极点。

例5 已知某离散系统的系统函数为H(z) =(2z+ 1)/( 3z^5 - 2z^4+ 1)试用MATLAB求出该系统的零极点,并画出零极点分布图,判断系统是否稳定。用roots()求得H(z)的零极点后，就可以用plot()函数绘制出系统的零极点图。

子函数

function ljdt(A,B)
p = roots(A);
%求系统极点
q= roots(B);
%求系统零点
p=p';
%将极点列向量转置为行向量
q=q';
%将零点列向量转置为行向量
x = max(abs([p q 1]));
%确定纵坐标范围
x=x+0.1;
y= x;
%确定横坐标范围
clf
hold on
axis([-x x -y y])
%确定坐标轴显示范围
w=0:pi/300:2*pi;
t= exp(i*w);
plot(t)            %画单位圆
axis( 'square')
plot([-x x],[0 0]) %画横坐标轴
plot([0 0],[- y y])%画纵坐标轴
text(0.1,x, 'JIm[z]')
text(y,1/10, 'Re[z]')
plot(real(p), imag(p), 'x')%画极点
plot(real(q), imag(q), 'o')
title('零极点图')%标注标题
hold off

主函数

a=[3 -2 0 0 0 1];
b=[2 1];
ljdt(a,b)
p = roots(a)
q = roots(b)
pa = abs(p)

例6 各种系统零极点图的实现

%绘制情况(a)系统零极点分布困及系统单位序列响应
z=0;      %定义系统零点位置
p= 0.25;  %定义系统极点位置
k=1;      %定义系统增益
subplot(221)
zplane(z, p)
grid on;
%绘制系统零极点分布图
subplot(222);
[ num, den]= zp2tf(z,p,k);%零极点模型转换为传递函数模型
impz( num, den)
%绘制系统单位序列响应时域波形
title( 'h(n)')
grid on; 
%定义标题
%绘制情况(b)系统零极点分布图及系统单位序列响应
p= 1;
subplot(223);
zplane(z, p)
grid on; 
[num, den] = zp2tf(z,p,k);
subplot(224);
impz(num, den)
title( 'b(n)')
grid on; 

4.离散系统差分函数求解

差分方程的解也分为通解与特解,通解是与方程初始状态有关的解,特解与外部输人有关,它描述系统在外部输人作用下的强迫运动。

例7 求解如下差分方程:
y(n)-0.5y(n-1)-0.45y(n-2) = 0.55x(n) +0.5x(n-1) - x(n- 2)，其中,x(n)=0.7^n u(n),初始状态y(-1)=1,y(- 2)=2,x(- 1)=2,x(-2)=3。

num=[0.55 0.5 -1];
den= [1  -0.5 -0.45];
x0=[2 3];y0=[1 2];
N= 50;
n=[0:N- 1]';
x=0.7.^n;
Zi= filtic(num, den, y0, x0);
[y,Zf]= filter(num, den, x,Zi);
plot(n,x,'r- ',n,y,'b--');
title('响应');
xlabel('n');ylabel('x(n)一y(n)'); 
legend('输入x','输出 y');
grid;

例8 编制程序求解下列两个系统差分方程的单位脉冲响应,并绘出其图形。
y[n]+0.75y[n-1]+0.125y[n-2]= x[n]-x[n- 1]

clc;
N= 32;
x_delta= zeros(1,N);
x_delta(1)=1;
p = [1, -1,0]
d=[1,0.75,0.125];
hl_delta= filter(p,d,x_delta);
subplot(211);
stem(0:N- 1,hl_delta, 'r');hold off;
xlabel('方程1的单位脉冲响应');
x_unit= ones(1,N);
h1_unit= filter(p,d,x_unit);
subplot(212);
stem(0:N- 1,h1_unit, 'r');hold off;
xlabel('方程1的阶跃响应');

我用的matlab 版本是2018a的，可能不同版本有些函数的用法稍微有些不同。