Biomedical Signal processing matlab Zhongguo Liu Biomedical Engineering
Biomedical Signal processing matlab 信号处理函数 Zhongguo Liu Biomedical Engineering School of Control Science and Engineering, Shandong University 2020/12/3 1 Zhongguo Liu_Biomedical Engineering_Shandong Univ.
MATLAB与信号处理直接有关的 具箱 ( Toolbox) Signal Processing (信号处理 具箱) Wavelet (小波 具箱) Image Processing(图象处理 具箱) Higher-Order Spectral Analysis (高阶谱分析 具箱) 3
与信号处理间接有关的 具箱: Control System (控制系统) Communication (通信) System Identification (系统辨识) Statistics (统计) Neural Network (神经网络) 4
例: z=peaks; 5 surf(z);
Z变换 u例③ .已知一离散信号的Z变换式为 , 求出它所对应的离散信号f(k). MATLAB程序如下: syms k z Fz=2*z/(2*z-1); %定义Z变换表达式 fk=iztrans(Fz, k) %求反Z变换 u运行结果如下:fk = (1/2)^k u例④:求序列的Z变换. syms n hn=sym('kronecker. Delta(n, 1) + kronecker. Delta(n, 2)+ kronecker. Delta(n, 3)') Hz=ztrans(hn) Hz=simplify(Hz) u运行结果如下:Fz= (z^2 + z + 1)/z^3 13
例1: 系统 用极零分析大致画出幅频响应及相频响应: 解: 转移函数: b=[1 -4 4]; a=[1]; [H, w]=freqz(b, a, 128, 'whole', 1) Hr=abs(H); Hphase=angle(H); Hphaseu=unwrap(Hphase); subplot(311), plot(Hr) subplot(312), plot(Hphase) subplot(313), plot(Hphaseu) 21
例2: b=[1]; a=[0, 1]; [H, w]=freqz(b, a, 256, 'whole', 1); Hr=abs(H); Hphase=angle(H); 相位的卷绕 (wrapping) 解卷绕 22 Hphase 1=unwrap(Hphase);
例: 给定系统 求:极-零图 解:b=[. 1836. 7344 1. 1016. 7374. 1836]/100 a =[1 -3. 0544 3. 8291 -2. 2925. 55075] 单位抽样响应 zplane(b, a); 频率响应 [h, t]=impz(b, a, 40); stem(t, h, '. '); grid on; [H, w]=freqz(b, a, 256, 'whole', 1); Hr=abs(H); 23
24 zplane(b, a); 极-零图
[h, t]=impz(b, a, 40); stem(t, h, '. '); grid on; 单位抽样响应 25
[H, w]=freqz(b, a, 256, 'whole', 1); Hr=abs(H); subplot(311), plot(Hr) subplot(312), plot(Hphase) subplot(313), plot(Hphaseu) Hphase=angle(H); 26 频率响应
例6. 7. 1(例6. 5. 1) 设计 IIR LP DF, uclear all; fp=100; fs=300; Fs=1000; rp=3; rs=20; uwp=2*pi*fp/Fs; ws=2*pi*fs/Fs; u. Fs=Fs/Fs; % let Fs=1 u% Firstly to finish frequency prewarping ; uwap=tan(wp/2); was=tan(ws/2); % u[n, wn]=buttord(wap, was, rp, rs, 's') % Note: 's'! u[z, p, k]=buttap(n); % u[bp, ap]=zp 2 tf(z, p, k) % u[bs, as]=lp 2 lp(bp, ap, wap) % u% Note: s=(2/Ts)(z-1)/(z+1); Ts=1, that is 2 Fs=1, Fs=0. 5; 32
例6. 7. 1(例6. 5. 1) 设计 IIR LP DF, uclear all; uwp=. 2*pi; ws=. 6*pi; Fs=1000; rp=3; rs=20; % u% Firstly to finish frequency prewarping; uwap=2*Fs*tan(wp/2); was=2*Fs*tan(ws/2); u[n, wn]=buttord(wap, was, rp, rs, 's'); % Note: 's'! u[z, p, k]=buttap(n); u[bp, ap]=zp 2 tf(z, p, k); u[bs, as]=lp 2 lp(bp, ap, wap) uw 1=[0: 499]*2*pi; uh 1=freqs(bs, as, w 1); u 33 [bz, az]=bilinear(bs, as, Fs) % Note: z=(2/ts)(z-
例6. 7. 1(例6. 5. 1) 设计 IIR LP DF, uclear all; uwp=. 2*pi; ws=. 6*pi; Fs=1000; urp=3; rs=20; u[n, wn]=buttord(wp/pi, ws/pi, rp, rs); u[bz, az]=butter(n, wp/pi) u[bz 1, az 1]=butter(n, wn) u[h, w]=freqz(bz, az, 128, Fs); u[h 1, w 1]=freqz(bz 1, az 1, 128, Fs); uplot(w, abs(h), w 1, abs(h 1), 'g. '); grid on; 34
9.cheb 2 ord. m; 10. ellipord. m; 11. cheb 2 ap. m; 12. ellipap. m; 13. besselap. m; 对应 Cheby-II、 椭圆 IIR 滤波器 14. cheby 2. m; 15. ellip. m; 16. besself. m 17.impinvar. m 用冲激响应不变法实现频率转换; 44
例7. 1. 1. 设计低通 DF FIR, 令截止频率 = 0. 25π, 取 M= 10, 20,40,观察其幅频响应的特点. clear all; N=10; b 1=fir 1(N, 0. 25, boxcar(N+1)); b 3=fir 1(2*N, 0. 25, boxcar(2*N+1)); b 4=fir 1(4*N, 0. 25, boxcar(4*N+1)); M=128; h 1=freqz(b 1, 1, M); h 3=freqz(b 3, 1, M); h 4=freqz(b 4, 1, M); f=0: 0. 5/M: 0. 5 -0. 5/M; plot(f, abs(h 1), f, abs(h 3), f, abs(h 4)); grid; axis([0 0. 5 0 1. 2]) 53
clear all; N=40; n=0: N; 例7. 1. 2: 理想差 b 1=fir 1(N, 0. 25, boxcar(N+1)); 分器及其设计 b 2=fir 1(N, 0. 25, hamming(N+1)); win=hamming(N+1); for n=1: N+1 if (n-1 -N/2)==0; b 1(n)=0; else M=128; b 1(n)=(-1)^(n-1 -N/2)/(n-1 -N/2); h 1=freqz(b 1, 1, M); end h 2=freqz(b 2, 1, M); end % h=freqz(b, 1, M); for n=1: N+1 f=0: 0. 5/M: 0. 5 -0. 5/M; if (n-1 -N/2)==0; hd=2*pi*f; b 2(n)=0; plot(f, abs(h 1), f, abs(h 2), f, hd, 'k-') else b 2(n)=win(n)*(-1)^(n-1 -N/2)/(n-1 -N/2); end 54
例7. 4. 2: 设计多阻带滤波器,抽样频率500 Hz, 在 50 Hz、 100 Hz 及150 Hz处陷波。 250 Hz clear all; f=[0. 14. 18. 22. 26. 34. 38. 42. 46. 54. 58. 62. 66 1]; A=[1 1 0 0 1 1]; weigh 1=[8 1 8 1 8]; b 1=remez(64, f, A, weigh 1); [h 1, w 1]=freqz(b 1, 1, 256, 1); h 1=abs(h 1); h 1=20*log 10(h 1); subplot(211); plot(w 1, h 1); grid; axis([0 0. 5 -60 10]) title('N=32, weight=[8 1 8 1 8]', 'Font. Size', 14, 'Color', 'r') 59
例3. 9. 1令x(n)为一正弦加白噪声信号,长度为 500, h(n) 是用fir 1. m文件设计出的一个低通FIR滤波器,长度为 11. 试用fftfilt实现长序列的卷积 clear; % 用叠接相加法,计算滤波器系数h和输入信号x的卷积 % 其中h为 10阶hanning窗,x是带有高斯白噪的正弦信号 h=fir 1(10, 0. 3, hanning(11)); % h: is the impulse response N=500; p=0. 05; f=1/16; % of a low-pass filter. u=randn(1, N)*sqrt(p); % u: white noise s=sin(2*pi*f*[0: N-1]); % s: sine signal x=u(1: N)+s; % x: a long sequence; y=fftfilt(h, x); % y=x*h subplot(211) plot(x); subplot(212) plot(y); 62
63
clear; % 生成滤波器系数h和混有高斯白噪的正弦信号x h=fir 1(10, 0. 3, hanning(11)); N=500; p=0. 05; f=1/16; u=randn(1, N)*sqrt(p); % s=sin(2*pi*f*[0: N-1]); x=u(1: N)+s; % 将x分为长度为L的小段 L=20; M=length(h); y=zeros(1, N+M-1); tempy=zeros(1, M+L-1); temp. X=zeros(1, L); for k=0: N/L-1 tempx(1: L)=x(k*L+1: (k+1)*L); tempy=conv(tempx, h); y=y+[zeros(1, k*L), tempy, zeros(1, N-(k+1)*L)]; end subplot(211); plot(x) 64 subplot(212); plot(y(1: N))
65
• 例 4. 7. 2 设 x(n)是 两 个 正 弦 信 号 及 白 噪 声 的 叠 加 , 进 行 频谱分析。 clear all; % 产生两个正弦加白噪声; N=256; f 1=. 1; f 2=. 2; fs=1; a 1=5; a 2=3; w=2*pi/fs; x=a 1*sin(w*f 1*(0: N-1))+a 2*sin(w*f 2*(0: N-1))+randn(1, N); % 应用FFT 求频谱; subplot(3, 1, 1); plot(x(1: N/4)); f=-0. 5: 1/N: 0. 5 -1/N; X=fft(x); y=ifft(X); subplot(3, 1, 2); plot(f, fftshift(abs(X))); subplot(3, 1, 3); plot(real(y(1: N/4))); 68
69
• 程序 clear all; % 构造三个不同频率的正弦信号的叠加作为试验信号 N=128; f 1=8; f 2=8. 22; f 3=9; fs=40; stepf=fs/N; n=0: N-1; t=2*pi*n/fs; n 1=0: stepf: fs/2 -stepf; x=sin(f 1*t)+sin(f 2*t)+sin(f 3*t); M=N; W=exp(-j*2*pi/M); % A=1时的czt变换 A=1; Y 1=czt(x, M, W, A); subplot(311) plot(n 1, abs(Y 1(1: N/2))); grid on; 71
% DTFT Y 2=abs(fft(x)); subplot(312) plot(n 1, abs(Y 2(1: N/2))); grid on; % 详细构造A后的czt M=60; f 0=7. 2; DELf=0. 05; A=exp(j*2*pi*f 0/fs); W=exp(-j*2*pi*DELf/fs); Y 3=czt(x, M, W, A); n 2=f 0: DELf: f 0+(M-1)*DELf; subplot(313); plot(n 2, abs(Y 3)); grid on; 72
3.deconv. m :实现系统的反卷积,其调用格式: [q, r]=deconv(y, x); 也用来实现多项式除法。 clear all; k=0: 1: 7; x=k+1; h=ones(1, 4); y=conv(h, x); % y = x * h; [q, r]=deconv(y, x); % 由 y, x 作反卷积,求出 h; [q 1, r 1]=deconv(y, h); % 由 y, h 作反卷积,求出 x; subplot(321); stem(h, '. b'); ylabel(' h(n)'); subplot(322); stem(x, '. b'); ylabel(' x(n)'); subplot(323); stem(y, '. b'); ylabel(' y(n)'); subplot(325); stem(q, '. r'); ylabel(' q(n)'); subplot(326); stem(q 1, '. r'); ylabel(' q 1(n)'); clear all; % 实现多项式除法q=(x^3+1)/(x+1) y=[1 0 0 1]; x=[1 1]; [q, r]=deconv(y, x); q 75
- Slides: 78