書上的matlab代碼敲一遍

---

1.卷積實現

clear; clc; x=[1,0,-1,1,0,1]; kx=-2:3; subplot(2,2,1);stem(kx,x); xlabel('kx');ylabel('x'); h=[1,0,2,-1,1]; kh=-2:2; subplot(2,2,2);stem(kh,h); xlabel('kh');ylabel('h'); y=conv(x,h); k=kx(1)+kh(1):kx(end)+kh(end); subplot(2,2,3); stem(k,y); xlabel('n');ylabel('y(n)');

k=kx(1)+kh(1):kx(end)+kh(end);
這里運用了

---

2.計算差分方程

clear; clc; N=41; a=[0.8 -0.44 0.36 0.02]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)]; k=0:1:N-1; y=filter(a,b,x); stem(k,y); xlabel('n');ylabel('fudu');

---

3.計算系統函數DTFT

$$H(z)=\frac{0.8?0.44z^{?1}+0.36z^{?2}+0.02z^{?3}}{1+0.7z^{?1}?0.45z^{?2}?0.3z^{?3}}$$

用到freqz函數(Frequency response of digital filter)

clear; clc; k=256; num=[0.8 -0.44 0.36 0.02]; den=[1 0.7 -0.45 -0.6]; w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1); plot(w/pi,real(h));grid title('實部'); xlabel('\omega/ \pi');ylabel('幅度'); subplot(2,2,2); plot(w/pi,imag(h));grid title('虛部'); xlabel('\omega/ \pi');ylabel('幅度'); subplot(2,2,3); plot(w/pi,abs(h));grid title('幅度譜'); xlabel('\omega/ \pi');ylabel('幅度'); subplot(2,2,4); plot(w/pi,angle(h));grid title('相位譜'); xlabel('\omega/ \pi');ylabel('相位');

---

4.畫幅頻曲線

畫出$H(z)=\frac{1}{1?\alpha z^{?1}}$的幅頻曲線，$\alpha$取正負0.9

clear; clc; b=[1];a1=[1 -0.9];a2=[1 0.9]; w=linspace(0,2*pi,512); h1=freqz(b,a1,w); h2=freqz(b,a2,w); plot(w/pi,abs(h1),w/pi,abs(h2),':'); xlabel('\omega/ \pi');ylabel('幅度');

---

5.畫系統函數的零極點分布

---

實驗一 離散時間的信號和系統

1.典型序列的實現

畫出單位抽樣序列、單位階躍序列和實指數序列$x(n)=0.6^n$的圖形

impseq.m

function[x,n]=impseq(n0,n1,n2)%Generate x(n)=delta(n-n0);n1<=n0<=n2;if((n0<n1)|(n0>n2)|(n1>n2))%變量范圍判斷error('參數必須滿足n1<=n0<=n2')endn=[n1:n2];%生成n1到n2的整數序列x=[(n-n0)==0];%判斷是否為0,在為零的位置產生一個脈沖

[x,n]=impseq(n0,n1,n2)
n0 脈沖產生位置
n1 序列起始值
n2 序列結束值
x單位采樣序列
n n1到n2的整數序列

測試代碼

n0=0;n1=-9;n2=9; [x,n] = impseq(n0,n1,n2); stem(n,x),title('\delta(n)'),xlabel('n (sample number)'),ylabel('\delta(n)');

（2）單位階躍序列

stepseq.m

function [x,n] = stepseq(n0,n1,n2) % Generates x(n) = u(n-n0); n1 <= n0 <=n 2; % [x,n] = stepseq(n0,n1,n2) if ((n0<n1)|(n0>n2)|(n1>n2))%變量范圍判斷error('Arguments must satisfy n1 <= n0 <= n2') end n = [n1:n2];%產生n1到n2的序列 x = [(n-n0) >= 0];%大于0就為1，單位階躍序列

[x,n] = stepseq(n0,n1,n2)
n0 脈沖產生位置
n1 序列起始值
n2 序列結束值
x單位階躍序列
n n1到n2的整數序列

測試代碼

n0=5;n1=-9;n2=9; [x,n] = stepseq(n0,n1,n2); stem(n,x),title('2u(n-5)'),xlabel('n (sample number)'),ylabel('2u(n-5)') axis([-10 10 -1 3]);

（3）實指數序列$x(n)=0.6^n$ , 0<=n<=10

n=[0:10]; x=(0.6).^n; stem(n,x); title('實指數序列'); xlabel('n (sample number)'); ylabel('0.6.^n');

2、序列的運算

(1)已知x1=[3,2 ,1, 1.5 ,-1, 0.5 ,4 ,-2], n1=[-2:5]；x2=[4 ,-2 ,1.5 ,2 ,1, 3.5 ,-1, -3.5],n2=[3:10]求兩信號的積。

seqmult.m

function [y,n]=seqmult(x1,n1,x2,n2) n=min(min(n1),min(n2)):max(max(n1),max(n2)); y1=zeros(1,length(n)); y2=y1; y1(find((n>=(min(n1)))&(n<=max(n1))==1))=x1; y2(find((n>=(min(n2)))&(n<=max(n2))==1))=x2; y=y1.*y2

測試

x1=[3,2 ,1, 1.5 ,-1, 0.5 ,4 ,-2];ns1=-2; x2=[4 ,-2 ,1.5 ,2 ,1, 3.5 ,-1, -3.5];ns2=3; nf1=5; nf2=10; n1=ns1: nf1; n2=ns2: nf2; n=min(ns1,ns2):max(nf1,nf2); y1=zeros(1,length(n)); y2=y1; y1(find((n>=ns1)&(n<=nf1)==1))=x1; y2(find((n>=ns2)&(n<=nf2)==1))=x2; ym=y1.*y2; ya=y1+y2; subplot(221);stem(n1,x1,'.');ylabel('x1(n)');grid; subplot(223);stem(n2,x2,'.');xlabel('n');ylabel('x2(n)'); grid; subplot(222);stem(n,ya,'.');ylabel('y1(n)+y2(n)');grid; subplot(224);stem(n,ym,'.');xlabel('n');ylabel('y1(n)*y2(n)'); grid;

(2)已知x=(3,2, 0.5, 1, 1.5 ,-1.5),nx=(-2,3)； h=[3,1,1,1,2,2.5];nh=(-1,4) 求兩信號的卷積和。
convwthn.m

function[y,ny]=convwthn(x,nx,h,nh) ny1=nx(1)+nh(1); ny2=nx(end)+nh(end); y=conv(x,h);ny=[ny1:ny2];

測試程序

x=[3,2, 0.5, 1, 1.5 ,-1.5]; nx=[-2,3]; h=[3,1,1,1,2,2.5]; nh=[-1,4]; [y,ny]=convwthn(x,nx,h,nh); stem(ny,y,'.') ;xlabel('n');ylabel('y(n)'); grid;

3. 求差分方程的沖激響應

b=1; a=[1,-0.7,-0.6,1]; x=[1,0,0,0,0,0]; y=filter(b,a,x)

實驗二 離散傅里葉變換（DFT）

1.離散傅立葉變換
求x(n)序列的32點離散傅里葉變換DFT，其中x(n)= cos(nπ/7)
dft.m

function[xk]=dft(xn,N) n=[0:1:N-1]; %n的行向量 k=[0:1:N-1]; %k的行向量 WN=exp(-j*2*pi/N); %WN因子 nk=n'*k; %產生一個含nk值的N乘N矩陣 WNnk=WN.^nk;%DFT矩陣 xk=xn*WNnk;%DFT系數的行向量

測試

N=32; n=[0:1:N-1]; k=[0:1:N-1]; xn=cos(pi*n/7); [xk]=dft(xn,N); subplot(2,1,1); stem(n,xn);xlabel('n');ylabel('xn'); subplot(2,1,2); stem(k,abs(xk)); xlabel('k');ylabel(' abs(xk)');

2.DFT的性質
求有限長序列x(n)=8(0.4)n，的圓周移位xN(n)=x((n+m)) 20R20(n)
得出當m=5,m=9時的圓周移位圖。
cirshift.m

function y=cirshift(x,m,N) if length(x)>Nerror('N必須>=x的長度') end x=[x zeros(1,N-length(x))];%將x補零到長度為N n=[0:1:N-1]; n=mod(n-m,N); y=x(n+1);%matlab下標從1開始

測試程序

n=[0:1:19]; x=8*(0.4).^n; y1=cirshift(x,-5,20); y2=cirshift(x,-9,20); subplot(2,1,1); stem(n,y1);xlabel('n');ylabel('y1'); title('m=5的圓周移位'); subplot(2,1,2); stem(n,y2);xlabel('n');ylabel('y2'); title('m=9的圓周移位');

3.利用DFT計算線性卷積
求有限長序列x1(n)=[5,4,3,2,1]與x2(n)= [1,2,3]的不同點的圓周卷積(N=5
和7)；并計算x1(n)和 x2(n)的線性卷積，根據實驗結果討論線性卷積和圓周卷積的關系。

circonvt.m

function y=circonvt(x1,x2,N) if length(x1)>Nerror('N必須>=x1的長度') end if length(x2)>Nerror('N必須>=x2的長度') end x1=[x1 zeros(1,N-length(x1))]; x2=[x2 zeros(1,N-length(x2))]; m=[0:1:N-1]; x2=x2(mod(-m,N)+1); H=zeros(N,N); for n=1:1:NH(n,:)=cirshift(x2,n-1,N); end y=x1*H';

測試程序

x1=[5,4,3,2,1]; x2=[1,2,3]; y1=circonvt(x1,x2,5); n=0:4; subplot(2,1,1); stem(n,y1);xlabel('n');ylabel('y1'); title('5點圓周卷積'); y2=circonvt(x1,x2,7); n=0:6; subplot(2,1,2); stem(n,y2);xlabel('n');ylabel('y2'); title('7點圓周卷積');

4.求序列的fft

fs=500; N=512; n=0:N-1; t=n/fs; x=sin(100*pi*t)+sin(300*pi*t)+rand(size(t)); y=fft(x,N); m=abs(y); f=n*fs/N; subplot(2,1,1),plot(t,x); grid on; subplot(2,1,2),plot(f(1:N/2),m(1:N/2)); grid on;

實驗二 IIR數字濾波器的設計

1、低通巴特沃思模擬濾波器的設計。
設計一個低通巴特沃思模擬濾波器,指標如下:
通帶截止頻率:fp=3400Hz，通帶最大衰減:Rp =3db。
阻帶截止頻率:fs=4000Hz，通帶最大衰減:As =40db。
基本參數
fp=3400Hz，fs=4000Hz，Rp =3db，As =40db

fp=3400;wp=2*pi*fp; %輸入濾波器的通帶截止頻率 fs=4000;ws=2*pi*fs ; %輸入濾波器的阻帶截止頻率 Rp=3;As=40; %輸入濾波器的通阻帶衰減指標 N=ceil(log10((10^(As/10)-1)/(10^(Rp/10)-1))/(2*log10(ws/wp)));%濾波器階次。 wc=wp/((10^(Rp/10)-1)^(1/(2*N))) ; %3dB截止頻率 [z0,p0,k0]=buttap(N); %字母后加0表明這是原型濾波器的各指標，而不是所求的濾波器的 p=p0*wc;a=real(poly(p)); %以下求非歸一化低通濾波器的分子、分母多項式系數 k=k0*wc^N;b0= real(poly(z0)); b=k*b0; [H,w]=freqs(b,a); plot(w/(2*pi)/1000,20*log10(abs(H))); xlabel('f(KHz)');ylabel('dB');axis([0,12,-50,1]) set(gca,'xtickmode','manual','xtick',[0,1,2,3,4,5,6]); set(gca,'ytickmode','manual','ytick',[-50,-40,-30,-20,-10,0]); grid;

首先根據基本參數確定濾波器階次N

N=ceil(log10((10^{(As/10)-1)/(10}(Rp/10)-1))/(2*log10(ws/wp)));%濾波器階次。
ceil函數是向右取整

然后計算通帶截止頻率歸一化參數

wc=wp/((10^(Rp/10)-1)(1/(2*N))) ; %3dB截止頻率

然后查表確定巴特沃斯低通原型濾波器

[z0,p0,k0]=buttap(N); %字母后加0表明這是原型濾波器的各指標，而不是所求的濾波器的
buttap()根據階數返回原型濾波器參數

P0是a參數，z0是空，k0是1

這三句暫時沒弄懂
p=p0wc;a=real(poly§); %以下求非歸一化低通濾波器的分子、分母多項式系數
k=k0wc^N;b0= real(poly(z0));
b=k*b0;

求出頻率函數
[H,w]=freqs(b,a);
畫出圖像

可以看到通帶截止頻率大概為為3400hz，阻帶截止頻率大概為4000hz

2、模擬低通轉換為數字低通濾波器
已知一模擬濾波器的系統函數為

分別用脈沖響應不變法和雙線性變換法將轉換為數字濾波器系統函數，并繪出和的幅頻響應曲線。分別取Fs=1000Hz、500Hz，觀察脈沖響應不變法中存在的頻率混疊失真和雙線性變換法中存在的非線性失真。

clear;close all b=1000;a=[1,1000]; w=[0:1000*2*pi]; %設定模擬頻率 [hf,w]=freqs(b,a,w); %計算模擬濾波器頻響函數 subplot(1,3,1) plot(w/2/pi,abs(hf));grid;title('(a)模濾幅頻') xlabel('f(Hz)');ylabel('幅度'); Fs0=[1000,500]; for m=1:2Fs=Fs0(m)[d,c]=impinvar(b,a,Fs) %用impinvar函數實現離散化[f,e]=bilinear(b,a,Fs) %用bilinear函數實現離散化wd=[0:512]*pi/512; %設定數字歸一化頻率hw1=freqz(d,c,wd); %計算數字濾波頻響函數hw2=freqz(f,e,wd);% 畫出數字濾波器幅頻特性subplot(1,3,2);plot(wd/pi,abs(hw1)/abs(hw1(1)));grid;title('(b)數濾幅頻(impinvar)');hold on subplot(1,3,3);plot(wd/pi,abs(hw2)/abs(hw2(1)));grid;title('(c)數濾幅頻(bilinear)');hold on; end

3.濾波器設計工具使用
命令行輸入filterbuilder

選擇lowpass

或者
輸入filter

按照要求設置參數

最后Design Filter生成濾波器

2.要求用SPTOOL工具箱選擇Butterworth IIR算法設計一個采樣頻率2000Hz，通邊界頻率：f1=100Hz，f 2=300Hz；通帶紋波指標：Rp=3dB；阻帶衰減：Rs=20dB的IIR低通濾波器。

改好參數，生成

完成
接下來用一個帶有噪聲的信號測試一下濾波器

首先生成濾波器代碼


fil.m

function Hd = fil %FIL Returns a discrete-time filter object.% MATLAB Code % Generated by MATLAB(R) 9.0 and the Signal Processing Toolbox 7.2. % Generated on: 09-Dec-2019 18:44:52% Butterworth Lowpass filter designed using FDESIGN.LOWPASS.% All frequency values are in Hz. Fs = 2000; % Sampling FrequencyFpass = 100; % Passband Frequency Fstop = 300; % Stopband Frequency Apass = 3; % Passband Ripple (dB) Astop = 20; % Stopband Attenuation (dB) match = 'stopband'; % Band to match exactly% Construct an FDESIGN object and call its BUTTER method. h = fdesign.lowpass(Fpass, Fstop, Apass, Astop, Fs); Hd = design(h, 'butter', 'MatchExactly', match);% [EOF]

首先，建一個50hz和400hz混合的信號

clear;clc; f=1000;%采樣率fs為2000，自定義信號最高為1000，這里設置為1000 N=1024;%采樣點1024，方便fft變換 n=0:N-1; t=n/f; x=sin(2*pi*50*t)+sin(2*pi*400*t); subplot(2,1,1); plot(t,x);xlabel('時域'); y=fft(x); subplot(2,1,2); plot(t(1:512),abs(y(1:512)));xlabel('頻域');

然后用上述濾波器進行濾波

Hd = fil; x2=filter(Hd,x);%自定義濾波器濾波 subplot(2,1,1); plot(t,x2);xlabel('時域'); y2=fft(x2); subplot(2,1,2); plot(t(1:512),abs(y2(1:512)));xlabel('頻域');

可以發現，高頻信號已經被濾除

全部代碼

clear;clc; f=1000;%采樣率fs為2000，自定義信號最高為1000，這里設置為1000 N=1024;%采樣點1024，方便fft變換 n=0:N-1; t=n/f; x=sin(2*pi*50*t)+sin(2*pi*400*t); subplot(2,2,1); plot(t,x);xlabel('時域'); y=fft(x); subplot(2,2,2); plot(t(1:512),abs(y(1:512)));xlabel('頻域');Hd = fil; x2=filter(Hd,x);%自定義濾波器濾波 subplot(2,2,3); plot(t,x2);xlabel('時域'); y2=fft(x2); subplot(2,2,4); plot(t(1:512),abs(y2(1:512)));xlabel('頻域');

實驗二 FIR數字濾波器的設計

1.用各種窗函數法設計FIR數字濾波器
分別用矩形窗和Hanming窗設計線性相位FIR低通濾波器。要求通帶截止頻率，單位脈沖響應h(n)的長度N=21。給出h(n)及其幅頻特性曲線，并比較不同窗函數時的的過渡帶寬度和阻帶最小衰減。
用窗函數法設計FIR數字濾波器時，先求出相應的理想濾波器（現為理想低通）單位脈沖響應，再根據阻帶最小衰減選擇合適的窗函數w(n)，最后得到FIR濾波器的單位脈沖響應。
現，N=21。所以線性相位理想低通濾波器的單位脈沖響應為：

信號處理工具箱中有窗生成函數boxcar，hamming，hanning，blackman等。wn=boxcar(m)產生長度為m的矩形窗函數列向量wn。其他窗函數的調用格式相同。

%用窗函數法設計FIR低通濾波器 N=21;wc=pi/4; %理想低通濾波器參數 n=0:N-1;r=(N-1)/2; hdn=sin(wc*(n-r))/pi./(n-r+eps); %計算理想低通單位脈沖響應 if rem(N,2)~=0hdn(r+1)=wc/pi; end %N為奇數時,處理n=r時的0/0型 wn1=boxcar(N); %矩形窗 hn1=hdn.*wn1'; %加窗 subplot(2,2,1) stem(n,hn1, '.') xlabel('n'),ylabel('hn1'); title('矩形窗') hw1=fft(hn1,512); w1=2*[0:511]/512; subplot(2,2,2) plot(w1,20*log10(abs(hw1))); grid wn2=hamming(N); %hamming窗 hn2=hdn.*wn2'; %加窗 subplot(2,2,3) stem(n,hn2, '.') xlabel('n'),ylabel('hn2'); title('海明窗') hw2=fft(hn2,512); w2=2*[0:511]/512; subplot(2,2,4) plot(w2,20*log10(abs(hw2))); grid

2.2. 用函數fir1和fir2設計FIR濾波器
MATLAB信號處理工具箱提供了基于窗函數法的FIR濾波器的設計函數fir1和fir2，利用它們可以使設計更簡單。
（1）fir1
功能：基于窗函數的FIR濾波器設計-----標準頻率響應形狀
格式：b=fir1(N,wc,’ftype’,window)
說明：標準頻率響應指所設計的預期特性為理想頻率響應：包括低通、帶通、高通或帶阻特性。ftype和window可以默認為低通和海明窗。如b=fir1(N,wc)可得到截止頻率為wc且滿足線性相位條件的N階FIR低通濾波器，window默認選用hamming窗。當wc=[wc1,wc2]時，得到的是通帶為的帶通濾波器。
b=fir1(N,wc,’ftype’)可設計高通和帶阻隔濾波器。
當ftype=high時，設計高通濾波器。
當ftype=stop時，設計帶阻濾波器。
注意：在設計高通和帶阻隔濾波器時，階數N只能取偶數。
window默認為為hamming窗。可用的其他窗函數有boxcar、hanning、bartlett、blackman窗等。
如：b=fir1(N,wc,bartlett(N+1))為使用bartlett窗進行設計。
（2）fir2
功能：基于窗函數的FIR濾波器設計-------任意頻率響應形狀
格式：b=fir2(N,f,m,window)
說明：fir2函數用于設計具有任意頻率響應形狀的加窗線性相位FIR數字濾波器，其幅頻特性由頻率點向量m給出。，要求f為單增向量，而且從0開始，以1結束，1表示數字頻率。M與f等長度，m(k)表示頻點f(k)的幅頻響應值。plot(f,m)命令可畫出期望逼近的幅頻響應曲線。N和window與fir1中的相同。
要求：改用fir1函數完成用矩形窗和Hamming窗設計線性相位FIR低通濾波器。要求歸一化截止頻率，單位脈沖響應h(n)的長度N=21。給出h(n)及其幅頻特性曲線，

1.矩形窗設計的線性相位FIR低通濾波器

n=21; %設濾器的階數為21 Wn=0.25;%歸一化截止頻率為0.25 window=rectwin(22);%矩形窗 b=fir1(n,Wn,window); freqz(b);

2.Hamming窗設計的線性相位FIR低通濾波器

n=21; %設濾器的階數為21 Wn=0.25;%歸一化截止頻率為0.25 window=hamming(22);%矩形窗 b=fir1(n,Wn,window); freqz(b);

總結

以上是生活随笔為你收集整理的DSP实验的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。