数字微分器的MATLAB设计

%hamming窗
M=21;alpha=(M-1)/2;n=0:M-1;
hd=(cos(pi*(n-alpha)))./(n-alpha);hd(alpha+1)=0;
w_ham=(hamming(M))'; h=hd.*w_ham; [Hr,w,P,L]=Hr_Type3(h);
%plots
subplot(221);stem(n,hd);
title('Ideal Impulse Response') axis([-1 M -1.2 1.2]);
xlabel('n');ylabel('hd(n)') subplot(222);stem(n,w_ham);
title('Hamming Window') axis([-1 M 0 1.2]);xlabel('n');ylabel('w(n)') 
subplot(223);stem(n,h);
title('Actual Impulse Response') axis([-1 M -1.2 1.2]);xlabel('n');ylabel('h(n)')
subplot(224);plot(w/pi,Hr/pi);
title('Amplitude Response');grid;
axis([0 1 0 1]);xlabel('frenquency in pi units');ylabel('slope in pi units');


子函数：function [Hr,w,c,L] = Hr_Type3(h);
% Computes Amplitude response of Type-3 LP FIR filter
% -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  ---
% [Hr,w,c,L] = Hr_Type3(h)

% Hr = Amplitude Response
%  w = frequencies between [0 pi] over which Hr is computed
%  c = Type-3 LP filter coefficients
%  L = Order of Hr
%  h = Type-3 LP impulse response
%
M = length(h);
L = (M-1)/2;
c = [2*h(L+1:-1:1)];
n=[0:1:L];
w = [0:1:500]'*pi/500;
Hr = sin(w*n)*c';

%blackman窗
M=21;alpha=(M-1)/2;n=0:M-1;
hd=(cos(pi*(n-alpha)))./(n-alpha);
hd(alpha+1)=0; w_black=(blackman(M))';
h=hd.*w_black; [Hr,w,P,L]=Hr_Type3(h);
%plots figure;
subplot(221);stem(n,hd);
title('Ideal Impulse Response') axis([-1 M -1.2 1.2]);
xlabel('n');ylabel('hd(n)') 
subplot(222);
stem(n,w_black);
title('Blackman Window') axis([-1 M 0 1.2]);xlabel('n');ylabel('w(n)') 
subplot(223);
stem(n,h);
title('Actual Impulse Response')
axis([-1 M -1.2 1.2]);xlabel('n');ylabel('h(n)') 
subplot(224);plot(w/pi,Hr/pi);
title('Amplitude Response');grid; 
axis([0 1 0 1]);xlabel('frenquency in pi units');ylabel('slope in pi units');

%频率采样法
M=33;alpha=(M-1)/2;
Dw=2*pi/M;l=0:M-1;wl=Dw*l; 
k1=0:floor((M-1)/2);k2=floor((M-1)/2)+1:M-1; 
Hrs=[j*Dw*k1,-j*Dw*(M-k2)];
angH=[-alpha*Dw*k1,alpha*Dw*(M-k2)];
H=Hrs.*exp(j*angH);h=real(ifft(H,M));
[Hr,ww,a,P]=Hr_Type3(h);
%plots figure;
subplot(211);stem(0:M-1,h);
title('Actual Impulse Response') axis([-1 M -0.5 1.5]);
xlabel('n');ylabel('h(n)') subplot(212);
plot(ww/pi,Hr/pi);title('Amplitude Response');grid;
axis([0 1 -0.5 1.5]);xlabel('frenquency in pi units');
ylabel('slope in pi units');