数字信号处理实验指导(6)

function [b,a] = afd_elip(Wp,Ws,Rp,As); % Analog Lowpass Filter Design: Elliptic % -------------------------------------- % [b,a] = afd_elip(Wp,Ws,Rp,As);

% b = Numerator coefficients of Ha(s) % a = Denominator coefficients of Ha(s)

% Wp = Passband edge frequency in rad/sec; Wp > 0

% Ws = Stopband edge frequency in rad/sec; Ws > Wp > 0 % Rp = Passband ripple in +dB; (Rp > 0)

% As = Stopband attenuation in +dB; (As > 0) %

if Wp <= 0

error('Passband edge must be larger than 0') end

if Ws <= Wp

error('Stopband edge must be larger than Passband edge') end

if (Rp <= 0) | (As < 0)

error('PB ripple and/or SB attenuation ust be larger than 0') end

ep = sqrt(10^(Rp/10)-1); A = 10^(As/20); OmegaC = Wp; k = Wp/Ws;

k1 = ep/sqrt(A*A-1);

capk = ellipke([k.^2 1-k.^2]); % Version 4.0 code

capk1 = ellipke([(k1 .^2) 1-(k1 .^2)]); % Version 4.0 code N = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1)));

fprintf('\\n*** Elliptic Filter Order = %2.0f \\n',N) [b,a]=u_elipap(N,Rp,As,OmegaC); function w_black = Blackman(M); % M-point Blackman window % ----------------------- % w_black = Blackman(M); %

M1 = M-1; m = [0:1:M1];

w_black = abs(0.42 - 0.5*cos(2*pi*m'/(M1)) + 0.08*cos(4*pi*m'/(M1)));

function [b,a] = cas2dir(b0,B,A); % CASCADE-to-DIRECT form conversion % ---------------------------------

26

% [b,a] = cas2dir(b0,B,A)

% b = numerator polynomial coefficients of DIRECT form % a = denominator polynomial coefficients of DIRECT form % b0 = gain coefficient

% B = K by 3 matrix of real coefficients containing bk's % A = K by 3 matrix of real coefficients containing ak's %

[K,L] = size(B); b = [1]; a = [1]; for i=1:1:K

b=conv(b,B(i,:)); a=conv(a,A(i,:)); end

b = b*b0;

function y = casfiltr(b0,B,A,x);

% CASCADE form realization of IIR and FIR filters % ----------------------------------------------- % y = casfiltr(b0,B,A,x); % y = output sequence

% b0 = gain coefficient of CASCADE form

% B = K by 3 matrix of real coefficients containing bk's % A = K by 3 matrix of real coefficients containing ak's % x = input sequence %

[K,L] = size(B); N = length(x); w = zeros(K+1,N); w(1,:) = x; for i = 1:1:K

w(i+1,:) = filter(B(i,:),A(i,:),w(i,:)); end

y = b0*w(K+1,:);

function [b,a] = cheb1hpf(wp,ws,Rp,As)

% IIR Highpass filter design using Chebyshev-1 prototype % [b,a] = cheb1hpf(wp,ws,Rp,As)

% b = Numerator polynomial of the highpass filter % a = Denominator polynomial of the highpass filter % wp = Passband frequency in radians % ws = Stopband frequency in radians % Rp = Passband ripple in dB

27

% As = Stopband attenuation in dB %

% Determine the digital lowpass cutoff frequecies: wplp = 0.2*pi;

alpha = -(cos((wplp+wp)/2))/(cos((wplp-wp)/2));

wslp = angle(-(exp(-j*ws)+alpha)/(1+alpha*exp(-j*ws))); %

% Compute Analog lowpass Prototype Specifications: T = 1; Fs = 1/T;

OmegaP = (2/T)*tan(wplp/2); OmegaS = (2/T)*tan(wslp/2);

% Design Analog Chebyshev Prototype Lowpass Filter: [cs,ds] = afd_chb1(OmegaP,OmegaS,Rp,As);

% Perform Bilinear transformation to obtain digital lowpass [blp,alp] = bilinear(cs,ds,Fs);

% Transform digital lowpass into highpass filter Nz = -[alpha,1]; Dz = [1,alpha]; [b,a] = zmapping(blp,alp,Nz,Dz); function [xec, xoc] = circevod(x)

% signal decomposition into circular-even and circular-odd parts% --------------------------------------------------------------% [xec, xoc] = circecod(x) %

if any(imag(x) ~= 0)

error('x is not a real sequence') end

N = length(x); n = 0:(N-1); xec = 0.5*(x + x(mod(-n,N)+1)); xoc = 0.5*(x - x(mod(-n,N)+1));

function y = circonvt(x1,x2,N)

% N-point circular convolution between x1 and x2: (time-domain) % ------------------------------------------------------------- % [y] = circonvt(x1,x2,N)

% y = output sequence containing the circular convolution % x1 = input sequence of length N1 <= N % x2 = input sequence of length N2 <= N % N = size of circular buffer

% Method: y(n) = sum (x1(m)*x2((n-m) mod N))

28

% Check for length of x1 if length(x1) > N

error('N must be >= the length of x1') end

% Check for length of x2 if length(x2) > N

error('N must be >= the length of 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:N

H(n,:) = cirshftt(x2,n-1,N); end

y = x1*H';

function y = cirshftt(x,m,N)

% Circular shift of m samples wrt size N in sequence x: (time domain) % ------------------------------------------------------------------- if length(x) > N

error('N must be >= the length of x') end

x = [x zeros(1,N-length(x))]; n = [0:1:N-1]; n = mod(n-m,N); y = x(n+1);

function [y,ny] = conv_m(x,nx,h,nh)

% Modified convolution routine for signal processing % -------------------------------------------------- % [y,ny] = conv_m(x,nx,h,nh) % y = convolution result % ny = support of y

% x = first signal on support nx % nx = support of x

% h = second signal on support nh

29

% nh = support of h %

nyb = nx(1)+nh(1); nye = nx(length(x)) + nh(length(h)); ny = [nyb:nye]; y = conv(x,h);

function [y,H]=conv_tp(h,x)

% Linear Convolution using Toeplitz Matrix % ---------------------------------------- % [y,H] = conv_tp(h,x)

% y = output sequence in column vector form

% H = Toeplitz matrix corresponding to sequence h so that y = Hx % h = Impulse response sequence in column vector form % x = input sequence in column vector form %

Nx = length(x); Nh = length(h); hc=[h; zeros(Nx-1, 1)]; hr=[h(1),zeros(1,Nx-1)]; H=toeplitz(hc,hr); y=H*x;

function I = cplxcomp(p1,p2) % I = cplxcomp(p1,p2)

% Compares two complex pairs which contain the same scalar elements % but (possibly) at differrent indices. This routine should be % used after CPLXPAIR routine for rearranging pole vector and its % corresponding residue vector. % p2 = cplxpair(p1) % I=[];

for j=1:1:length(p2) for i=1:1:length(p1)

if (abs(p1(i)-p2(j)) < 0.0001) I=[I,i]; end end end I=I';

function [d1,d2] = db2delta(Rp,As)

% Conversion from Relative dB specs to Absolute delta specs. % ---------------------------------------------------------- % [d1,d2] = db2delta(Rp,As)

30

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库数字信号处理实验指导(6)在线全文阅读。