magXw=abs(Xw);angXw=angle(Xw); magXk=abs(Xk);angXk=angle(Xk);
subplot(2,2,1);plot(w/pi,magXw);grid
xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,3);plot(w/pi,angXw);grid
xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度') subplot(2,2,2);stem(k,magXk);grid
xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,4);stem(k,angXk);grid
xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度')
function[X]=dtft(x,n,w) X=x*(exp(-j).^(n'*w));
function[y]=dtftplot(w,X) function[y]=dtftplot(w,X) magX=abs(X);angX=angle(X); realX=real(X);imagX=imag(X);
subplot(2,2,1);plot(w/pi,magX);grid
xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,3);plot(w/pi,angX);grid
xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度') subplot(2,2,2);plot(w/pi,realX);grid
xlabel('以pi为单位的频率');title('实部部分');ylabel('实部') subplot(2,2,4);plot(w/pi,imagX);grid
xlabel('以pi为单位的频率');title('虚部部分');ylabel('虚部')
function [xe, xo, m] = evenodd(x,n)
% Real signal decomposition into even and odd parts % ------------------------------------------------- % [xe, xo, m] = evenodd(x,n) %
if any(imag(x) ~= 0)
error('x is not a real sequence') end
m = -fliplr(n);
m1 = min([m,n]); m2 = max([m,n]); m = m1:m2; nm = n(1)-m(1); n1 = 1:length(n); x1 = zeros(1,length(m)); x1(n1+nm) = x; x = x1; xe = 0.5*(x + fliplr(x)); xo = 0.5*(x - fliplr(x));
36
function [db,mag,pha,w] = freqs_m(b,a,wmax);
% Computation of s-domain frequency response: Modified version % ------------------------------------------------------------ % [db,mag,pha,w] = freqs_m(b,a,wmax);
% db = Relative magnitude in db over [0 to wmax] % mag = Absolute magnitude over [0 to wmax]
% pha = Phase response in radians over [0 to wmax]
% w = array of 500 frequency samples between [0 to wmax] % b = Numerator polynomial coefficents of Ha(s) % a = Denominator polynomial coefficents of Ha(s)
% wmax = Maximum frequency in rad/sec over which response is desired %
w = [0:1:500]*wmax/500; H = freqs(b,a,w); mag = abs(H);
db = 20*log10((mag+eps)/max(mag)); pha = angle(H);
function [db,mag,pha,grd,w] = freqz_m(b,a); % Modified version of freqz subroutine % ------------------------------------ % [db,mag,pha,grd,w] = freqz_m(b,a);
% db = Relative magnitude in dB computed over 0 to pi radians % mag = absolute magnitude computed over 0 to pi radians % pha = Phase response in radians over 0 to pi radians % grd = Group delay over 0 to pi radians
% w = 501 frequency samples between 0 to pi radians % b = numerator polynomial of H(z) (for FIR: b=h) % a = denominator polynomial of H(z) (for FIR: a=[1]) %
[H,w] = freqz(b,a,1000,'whole');
H = (H(1:1:501))'; w = (w(1:1:501))'; mag = abs(H);
db = 20*log10((mag+eps)/max(mag)); pha = angle(H);
% pha = unwrap(angle(H)); grd = grpdelay(b,a,w); % grd = diff(pha); % grd = [grd(1) grd];
% grd = [0 grd(1:1:500); grd; grd(2:1:501) 0]; % grd = median(grd)*500/pi; function [xn] = idfs(Xk,N)
37
% Computes Inverse Discrete Fourier Series % ---------------------------------------- % [xn] = idfs(Xk,N)
% xn = One period of periodic signal over 0 <= n <= N-1 % Xk = DFS coeff. array over 0 <= k <= N-1 % N = Fundamental period of Xk %
n = [0:1:N-1]; % row vector for n k = [0:1:N-1]; % row vecor for k WN = exp(-j*2*pi/N); % Wn factor
nk = n'*k; % creates a N by N matrix of nk values
WNnk = WN .^ (-nk); % IDFS matrix
xn = (Xk * WNnk)/N; % row vector for IDFS values
function [b,a] = u_buttap(N,Omegac);
% Unnormalized Butterworth Analog Lowpass Filter Prototype % -------------------------------------------------------- % [b,a] = u_buttap(N,Omegac);
% b = numerator polynomial coefficients of Ha(s) % a = denominator polynomial coefficients of Ha(s) % N = Order of the Butterworth Filter % Omegac = Cutoff frequency in radians/sec %
[z,p,k] = buttap(N); p = p*Omegac; k = k*Omegac^N; B = real(poly(z)); b0 = k; b = k*B;
a = real(poly(p));
function [b,a] = u_chb1ap(N,Rp,Omegac);
% Unnormalized Chebyshev-1 Analog Lowpass Filter Prototype % -------------------------------------------------------- % [b,a] = u_chb1ap(N,Rp,Omegac);
% b = numerator polynomial coefficients % a = denominator polynomial coefficients % N = Order of the Elliptic Filter % Rp = Passband Ripple in dB; Rp > 0 % Omegac = Cutoff frequency in radians/sec %
[z,p,k] = cheb1ap(N,Rp); a = real(poly(p));
38
aNn = a(N+1); p = p*Omegac;
a = real(poly(p)); aNu = a(N+1); k = k*aNu/aNn; b0 = k;
B = real(poly(z)); b = k*B;
function [b,a] = u_chb2ap(N,As,Omegac);
% Unnormalized Chebyshev-2 Analog Lowpass Filter Prototype% --------------------------------------------------------% [b,a] = u_chb2ap(N,As,Omegac);
% b = numerator polynomial coefficients % a = denominator polynomial coefficients % N = Order of the Elliptic Filter % As = Stopband Ripple in dB; As > 0 % Omegac = Cutoff frequency in radians/sec %
[z,p,k] = cheb2ap(N,As); a = real(poly(p)); aNn = a(N+1); p = p*Omegac;
a = real(poly(p)); aNu = a(N+1);
b = real(poly(z)); M = length(b); bNn = b(M); z = z*Omegac;
b = real(poly(z)); bNu = b(M);
k = k*(aNu*bNn)/(aNn*bNu); b0 = k; b = k*b;
function [b,a] = u_elipap(N,Rp,As,Omegac);
% Unnormalized Elliptic Analog Lowpass Filter Prototype % ----------------------------------------------------- % [b,a] = u_elipap(N,Rp,As,Omegac);
% b = numerator polynomial coefficients % a = denominator polynomial coefficients % N = Order of the Elliptic Filter % Rp = Passband Ripple in dB; Rp > 0
% As = Stopband Attebuation in dB; As > 0
39
% Omegac = Cutoff frequency in radians/sec %
[z,p,k] = ellipap(N,Rp,As); a = real(poly(p)); aNn = a(N+1); p = p*Omegac;
a = real(poly(p)); aNu = a(N+1);
b = real(poly(z)); M = length(b); bNn = b(M); z = z*Omegac;
b = real(poly(z)); bNu = b(M);
k = k*(aNu*bNn)/(aNn*bNu); b0 = k; b = k*b;
function [bz,az] = zmapping(bZ,aZ,Nz,Dz)
% Frequency band Transformation from Z-domain to z-domain% -------------------------------------------------------% [bz,az] = zmapping(bZ,aZ,Nz,Dz) % performs:
% b(z) b(Z)|
% ---- = ----| N(z) % a(z) a(Z)|@Z = ---- % D(z) %
bzord = (length(bZ)-1)*(length(Nz)-1); azord = (length(aZ)-1)*(length(Dz)-1);
bz = zeros(1,bzord+1); for k = 0:bzord pln = [1]; for l = 0:k-1
pln = conv(pln,Nz); end
pld = [1];
for l = 0:bzord-k-1 pld = conv(pld,Dz); end
bz = bz+bZ(k+1)*conv(pln,pld); end
40