数字信号处理MATLAB笔记

94
0
2020年11月10日 09时03分

全文框架

 

1

 

1.函数笔记

 

2

3

 

fft():快速傅里叶变换
Y = fft(X,n,dim)
计算x的n点DFT,x长度不够n时补0,比n长时截短,dim为对x做傅里叶变换的维数。

 

ifft():傅里叶反变换
y = ifft(X,n)
计算X的n点DFT反变换

 

Stem():火柴杆图
stem(X,Y)
显示离散序列,X与Y大小必须相同,一一对应。如Y是一个长度为6的序列,x可以为0:5,若X省略,程序会自动把横坐标变为Y的大小。
若一个连续序列用stem显示,会变成实心的,因为离散序列每个点下都有一条竖线,因此最好用plot显示连续序列。

 

Abs():绝对值或幅值
Y = abs(X)
求绝对值或复数的幅度,可用于求幅度谱。

 

Angle():相位
P = angle(Z)
求复数Z的角度和相位,返回值在±π之间。

 

Conj():求共轭
ZC = conj(Z)

矩形信号:
矩形序列R(N):可以用zeros和ones构造,如:R20=[zeros(1,40),ones(1,20),zeros(1,41)];
连续矩形信号可用rectpuls()函数,y = rectpuls(t,w),t为时间序列,w为宽度,宽度以t=0为对称
单位脉冲信号:
只需在原点处设一个脉冲即可

 

Tripuls():三角波信号
y = tripuls(T,width,skew)
T为时间序列,width为三角波宽度,skew为斜度,范围0-1,默认为0,即对称三角波。
注意三角波以T=0对称,若要完整显示,时间序列要从负到正(如:T=-30:1:30)。
若要移动三角波,对T进行加减(如T+50表示向右移50单位)

 

Impulse():单位脉冲响应
Y=Impulse(sys,t)
Sys为系统函数,t为时间序列(Ts:dt:Te),y为单位阶跃响应,t省略则自动生成时间间隔

 

Step():单位阶跃响应
Y=step(sys,t)
Sys为系统函数,t为时间序列(Ts:dt:Te),y为单位阶跃响应,t省略则自动生成时间间隔
[y,t] = step(sys)
返回单位阶跃响应y,时间间隔t

 

Filter():一维滤波
y = filter(b,a,x)
a,b为滤波器的传递函数的分子分母,x为向量或矩阵,y为输出,大小与x相同。
也可以用来求系统响应,比如系统差分方程: ,则b和a为: ,x即为x(n),输出y为系统响应

 

Impz():求离散LTI系统的单位冲激响应
h=impz(b,a,k)
b为x的向量,a为y的向量,k为输出序列的取值范围。h为单位脉冲响应。

 

tf():系统函数
sys = tf(Numerator,Denominator)
创建系统函数,numerator为分子向量,denominator为分母向量
后面可加参数,如:‘Variable’,‘p’,意思是把p作为变量

 

Freqz():离散系统频率响应
[h,w] = freqz(b,a,n)
B为分子,a为分母,返回n个频率点的频率响应,放在h向量中,n个频率放在向量w中,这n个频率均匀设置在(0,pi)上,缺省时默认512

 

Freqs():连续系统频率响应
h = freqs(b,a,w)
B,a为分子分母,w为角频率向量,返回频率响应h

 

Ellipord():椭圆滤波器参数
[n,Wp] = ellipord(Wp,Ws,Rp,Rs,’s’)
Rp、Rs分别为通带最大波纹和阻带最小衰减;wp、ws分别为为通带边界频率和阻带边界频率,单位为rad/s。
‘s’表示模拟滤波器,缺省时该函数适用于数字滤波器,此时wp及ws均为0~1之间的数值,即频率与抽样频率的比值。返回最小阶数n和截止频率Wp。

 

编写N点循环卷积程序:
X_K1=fft(xn1,N);
X_K2=fft(xn2,N);
X_K3=X_K1.*X_K2;
x_n3=ifft(X_K3,N);
因fft可以自动补0,故可以利用其将两序列变成长度为N的DFT序列,再对应相乘,相当于循环卷积的效果。此程序可用来比较线性卷积(conv函数)和循环卷积的区别。

 

其他笔记:

 

4

 

5\\

2.1系统响应及系统稳定性

一、实验目的和步骤
目的:
(1)掌握 求系统响应的方法。
(2)掌握时域离散系统的时域特性。
(3)分析、观察及检验系统的稳定性。
步骤:

 

7

 

二、代码

 

close all
clear all 
clc

x1_n=[ones(1,8),zeros(1,90)];
x2_n=ones(1,100);
t=0:1:99;
h_n=impz([0.05 0.05],[1 -0.9],t);%系统单位脉冲响应
y1_n=conv(x1_n,h_n);
y2_n=conv(x2_n,h_n);
figure
subplot(311);stem(h_n);title('系统单位脉冲响应');
subplot(312);stem(y1_n(1:100));title('R8(n)系统响应');
subplot(313);stem(y2_n(1:100));title('u(n)系统响应');

h1_n=[ones(1,10)];
h2_n=[1 2.5 2.5 1];
y3_n=conv(x1_n,h1_n);
y4_n=conv(x1_n,h2_n);
figure
subplot(211);stem(y3_n(1:100));title('对h1(n)输出响应');
subplot(212);stem(y4_n(1:100));title('对h2(n)输出响应');


xsin=sin(0.014*t)+sin(0.4*t); %产生正弦信号
A=[1,-1.8237,0.9801];B=[1/100.49,0,-1/100.49];  %系统差分方程系数向量B和A
y5_n=filter(B,A,x2_n);    %谐振器对u(n)的响应
y6_n=filter(B,A,xsin);    %谐振器对x(n)的响应
figure
subplot(2,1,1);stem(y5_n);title('谐振器对u(n)的响应');
subplot(2,1,2);stem(y6_n);title('谐振器对正弦信号的响应');

 

三、运行结果

 

8

9

10

 

2.2 时域采样与频域采样

一、实验目的和步骤
目的:
验证时域采样和频域采样定理
步骤:

 

11

12

 

二、代码

 

close all
clear all 
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%5时域采样定理验证
Tp=64/1000;%观察时间
Fs1=1000;%采样频率,分别为1000,300,200Hz
Fs2=300;
Fs3=200;
T1=1/Fs1;%周期
T2=1/Fs2;
T3=1/Fs3;
M=64;
A=444.128;
a=pi*50*2^0.5;
w=pi*50*2^0.5;
n1=0:Tp*Fs1;
n2=0:Tp*Fs2;
n3=0:Tp*Fs3;
xa_t1=A*exp(-a*n1*T1).*sin(w*n1*T1);
xa_t2=A*exp(-a*n2*T2).*sin(w*n2*T2);
xa_t3=A*exp(-a*n3*T3).*sin(w*n3*T3);
Xa_k1=fft(xa_t1,M);
Xa_k2=fft(xa_t2,M);
Xa_k3=fft(xa_t3,M);
figure
subplot(321);stem(xa_t1);title('采样频率为1000Hz的时域离散信号');
subplot(322);stem(abs(Xa_k1));title('64点DFT');
subplot(323);stem(xa_t2);title('采样频率为300Hz的时域离散信号');
subplot(324);stem(abs(Xa_k2));title('64点DFT');
subplot(325);stem(xa_t3);title('采样频率为300Hz的时域离散信号');
subplot(326);stem(abs(Xa_k3));title('64点DFT');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%频域采样定理的验证

for i=0:13
    x_n(i+1)=i+1;
end
for i=14:26
    x_n(i)=27-i;
end
X_k32=fft(x_n,32);
X_k16=X_k32(1:2:32);
x_n32=ifft(X_k32,32);
x_n16=ifft(X_k16,16);

figure
subplot(321);plot(abs(fft(x_n)));title('原频谱');
subplot(322);stem(x_n);title('原序列');
subplot(323);stem(abs(X_k32));title('32点采样频谱');
subplot(324);stem(x_n32);title('32点IDFT');
subplot(325);stem(abs(X_k16));title('16点采样频谱');
subplot(326);stem(x_n16);title('32点采样频谱');

 

三、运行结果

 

13

14

 

2.3 用FFT对信号作频谱分析

一、实验目的和步骤
步骤:

 

15

 

二、代码

 

close all 
clear all
clc

x1_n=[ones(1,4)];
for i=1:4
    x2_n(i)=i;
end
for i=5:8
    x2_n(i)=9-i;
end

for i=1:4
    x3_n(i)=5-i;
end
for i=5:8
    x3_n(i)=i-2;
end

X1_k8=fft(x1_n,8);
X1_k16=fft(x1_n,16);
X2_k8=fft(x2_n,8);
X2_k16=fft(x2_n,16);
X3_k8=fft(x3_n,8);
X3_k16=fft(x3_n,16);

figure
subplot(321);stem(abs(X1_k8));title('x1(n)的8点DFT');
subplot(322);stem(abs(X1_k16));title('x1(n)的16点DFT');
subplot(323);stem(abs(X2_k8));title('x2(n)的8点DFT');
subplot(324);stem(abs(X2_k16));title('x2(n)的16点DFT');
subplot(325);stem(abs(X3_k8));title('x3(n)的8点DFT');
subplot(326);stem(abs(X3_k16));title('x3(n)的16点DFT');

n=1:32;
x4_n=cos(pi/4*n);
x5_n=cos(pi*n/4)+cos(pi*n/8);
X4_k8=fft(x4_n,8);
X4_k16=fft(x4_n,16);
X5_k8=fft(x5_n,8);
X5_k16=fft(x5_n,16);

figure
subplot(221);stem(abs(X4_k8));title('x4(n)的8点DFT');
subplot(222);stem(abs(X4_k16));title('x4(n)的16点DFT');
subplot(223);stem(abs(X5_k8));title('x5(n)的8点DFT');
subplot(224);stem(abs(X5_k16));title('x5(n)的16点DFT');

N1=16;
N2=32;
N3=64;
n1=1:1:N1;
n2=1:1:N2;
n3=1:1:N3;
Fs=64;
t1=n1./Fs;
t2=n2./Fs;
t3=n3./Fs;
x8_t1=cos(8*pi*t1)+cos(16*pi*t1)+cos(20*pi*t1);
X8_k1=fft(x8_t1,N1);
x8_t2=cos(8*pi*t2)+cos(16*pi*t2)+cos(20*pi*t2);
X8_k2=fft(x8_t2,N2);
x8_t3=cos(8*pi*t3)+cos(16*pi*t3)+cos(20*pi*t3);
X8_k3=fft(x8_t3,N3);
figure
subplot(311);stem(abs(X8_k1));title('x8(t)的16点DFT');
subplot(312);stem(abs(X8_k2));title('x8(t)的32点DFT');
subplot(313);stem(abs(X8_k3));title('x8(t)的64点DFT');

 

三、运行结果

 

16

17

18

 

2.4 IIR数字滤波器设计

一、实验目的和步骤
目的:
(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法;
(2)学会调用MATLAB信号处理工具箱中滤波器设计函数(或滤波器设计分析工具fdatool)设计各种IIR数字滤波器,学会根据滤波需求确定滤波器指标参数。
(3)掌握IIR数字滤波器的MATLAB实现方法。
(3)通过观察滤波器输入输出信号的时域波形及其频谱,建立数字滤波的概念。
步骤:

 

20

111

 

二、代码

 

close all
clear all
clc

%%%%%%%%%%%%%滤波器设计
st=mstg;   %调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号s
Fs=10000; T=1/Fs;  %采样频率
fp=280; fs=450;wp=2*fp/Fs;ws=2*fs/Fs; %低通(关于pi归一化)
   name='y_1(t)';LHL_filter(wp,ws,name,st,T,'low');
fp_L=440; fp_R=560;fs_L=275; fs_R=900;%带通
   wp=[2*fp_L/Fs,2*fp_R/Fs];ws=[2*fs_L/Fs,2*fs_R/Fs];
   name='y_2(t)'; LHL_filter(wp,ws,name,st,T,'bandpass');
fp=890; fs=600;wp=2*fp/Fs;ws=2*fs/Fs; %高通
   name='y_3(t)';LHL_filter(wp,ws,name,st,T,'high');

function LHL_filter(wp,ws,name,st,T,flag)
rp=0.1;rs=60; %DF指标(低通滤波器的通、阻带边界频)
[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp?
[B,A]=ellip(N,rp,rs,wp,flag); %调用ellip计算椭圆带通DF系统函数系数向量B和A
yt=filter(B,A,st); %滤波器软件实现
 figure; freqz(B,A);%B,A分别为系统函数分子,分母多项式系数向量?
 figure;tplot(st,yt,T,name);
end
function tplot(st,xn,T,name)
%时域序列连续曲线绘图函数
%?xn:信号数据序列,%?T为采样间隔
N=length(xn);n=0:N-1; t=n*T; Tp=N*T; f=n/Tp;
subplot(311);plot(t,xn);xlabel('t/s');ylabel(name);
axis([0,t(end),min(xn),1.2*max(xn)]);
subplot(312);stem(f,abs(fftshift(fft(st,N))));
title('原 st(t)的频谱');xlabel('f/Hz');ylabel('幅度');
subplot(313);stem(f,abs(fftshift(fft(xn,N))));
title('滤波后xn(t)的频谱');xlabel('f/Hz');ylabel('幅度');
end
function st=mstg
%产生信号序列向量st,并显示st的时域波形和频谱
%st=mstg 返回三路调幅信号相加形成的混合信号,长度N=1600
N=1600;Fs=10000;T=1/Fs;Tp=N*T; % N为信号st的长度。采样频率Fs=10kHz,Tp为采样时间
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路调幅信号的载波频率fc1=1000Hz,调制信号频率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路调幅信号的载波频率fc2=500Hz,调制信号频率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路调幅信号的载波频率fc3=250Hz,调制信号频率fm3=25Hz                  
xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第1路调幅信号
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第2路调幅信号
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第3路调幅信号
st=xt1+xt2+xt3;         %三路调幅信号相加
fxt=fft(st,N);          %计算信号st的频谱
subplot(211);plot(t,st);xlabel('t/s');ylabel('s(t)');
axis([0,Tp/8,min(st),max(st)]);title('(a) s(t)的波形');
subplot(212);stem(f,abs(fxt)./max(abs(fxt)));
axis([0,Fs/5,0.1,2]);
title('(频谱');xlabel('f/Hz');ylabel('幅度');
end

 

三、运行结果

 

21

22

23

24

25

26

112

 

113

 

2.5 FIR数字滤波器设计

一、实验目的和步骤
目的:
(1)掌握用窗函数法设计FIR数字滤波器的原理和方法。
(2)掌握用等波纹最佳逼近法设计FIR数字滤波器的原理和方法。
(3)掌握FIR滤波器的快速卷积实现原理。
(4)学会调用MATLAB函数设计与实现FIR滤波器。
步骤:

 

28

27

 

 

二、代码

 

clc
close all
clear all

xt=xtg;
Fs=1000;T=1/Fs;
fp=120;fs=150;
wp=2*fp*T;ws=2*fs*T;
%%%%%%%%%%%%%%%%%%%%低通滤波器
rp=0.1;rs=60; %DF指标
[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B,A]=ellip(N,rp,rs,wp,'low'); %调用ellip计算椭圆带通DF系统函数系数向量B和A
low_filter=freqz(B,A);
yt1=filter(B,A,xt); 
figure
subplot(311);plot(abs(low_filter));title('椭圆低通滤波器幅度谱');xlabel('f/Hz');axis([0,600,0,1.2]);
subplot(312);plot(angle(low_filter));title('椭圆低通滤波器相位谱');xlabel('f/Hz');
subplot(313);plot(20*log(abs(low_filter)));title('椭圆低通滤波器损耗函数');xlabel('f/Hz');
figure
subplot(311);plot(yt1);title('椭圆低通滤波后时域波形');axis([0,600,-1.2,1.2]);
subplot(312);plot(abs(fft(yt1)));title('幅频特性');
subplot(313);plot(angle(fft(yt1)));title('相频特性');

%%%%%%%%%%%%%%%%%窗函数法,逼近的函数为理想低通滤波器,使用凯塞窗
Bt=ws-wp;%过渡带宽度
alph=0.112*(rs-8.7);%凯塞窗参数a
M=ceil((rs-8)/2.285/Bt)%凯塞窗阶数
wc=(ws+wp)/2;%已对pi归一化
hn=fir1(M,wc,kaiser(M+1,alph));
[HK,w]=freqz(hn,1);
yt2=fftfilt(hn,xt,1000);
figure
subplot(411);plot(hn);title('凯塞窗h(n)波形');xlabel('n');
subplot(412);plot(w/pi,abs(HK));title('幅度谱');xlabel('w/\pi');
subplot(413);plot(w/pi,angle(HK));title('相位谱');xlabel('w/\pi');
subplot(414);plot(w/pi,20*log(abs(HK)));title('损耗函数');xlabel('w/\pi');

figure
subplot(311);plot(yt2);title('凯塞窗滤波后时域波形');axis([0,600,-1.2,1.2]);
subplot(312);plot(abs(fft(yt2)));title('幅频特性');
subplot(313);plot(angle(fft(yt2)));title('相频特性');

%%%%%%%%%%%%%%%%%%%%%%%%%等波纹法
f=[fp,fs];m=[1,0];
dat1=(10^(rp/20)-1)/(10^(rp/20)+1);dat2=10^(-rs/20);
rip=[dat1,dat2];
[M1,fo,mo,w]=remezord(f,m,rip,Fs);
 M1=M1+1
hn1=remez(M1,fo,mo,w);
[HK1,w]=freqz(hn1,1);
yt3=fftfilt(hn1,xt,1000);
figure
subplot(411);plot(hn);title('等波纹法滤波器h(n)波形');xlabel('n');
subplot(412);plot(w/pi,abs(HK1));title('幅度谱');xlabel('w/\pi');
subplot(413);plot(w/pi,angle(HK1));title('相位谱');xlabel('w/\pi');
subplot(414);plot(w/pi,20*log(abs(HK1)));title('损耗函数');xlabel('w/\pi');

figure
subplot(311);plot(yt3);title('等波纹法滤波后时域波形');axis([0,600,-1.2,1.2]);
subplot(312);plot(abs(fft(yt3)));title('幅频特性');
subplot(313);plot(angle(fft(yt3)));title('相频特性');

function xt=xtg
N=1000;Fs=1000;T=1/Fs;Tp=N*T;
t=0:T:(N-1)*T;
fc=Fs/10;
f0=fc/10;
mt=cos(2*pi*f0*t);
ct=cos(2*pi*fc*t);
xt=mt.*ct;
nt=2*rand(1,N)-1;
fp=150;fs=200;Rp=0.1;As=70;
fb=[fp,fs];m=[0,1];
dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];
[n,fo,mo,W]=remezord(fb,m,dev,Fs);
hn=remez(n,fo,mo,W);
yt=filter(hn,1,10*nt);

xt=xt+yt;
fst=fft(xt,N);k=0:N-1;f=k/Tp;
subplot(311);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)');
axis([0,Tp/5,min(xt),max(xt)]);title('(a)加噪声波形');
subplot(312);plot(f,abs(fst)/max(abs(fst)));grid;title('(b)信号加噪声频谱');
axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度');
end

 

三、运行结果

 

31

32

33

34

35

36

35

总结

通过matlab编程,更深刻掌握了数字信号处理的知识,very interesting

 

发表评论

后才能评论