数字信号处理(基于计算机的方法 第三版)matlab程序.doc

% Program 2_1 % Generation of the ensemble average % R = 50; m = 0:R-1; s = 2*m.*(0.9.m); % Generate the uncorrupted signal d = rand(R,1)-0.5; % Generate the random noise x1 = s+d; stem(m,d); xlabel(Time index n);ylabel(Amplitude); title(Noise); pause for n = 1:50; d = rand(R,1)-0.5; x = s + d; x1 = x1 + x; end x1 = x1/50; stem(m,x1); xlabel(Time index n);ylabel(Amplitude); title(Ensemble average);% Program 2_2 % Generation of complex exponential sequence % a = input(Type in real exponent = ); b = input(Type in imaginary exponent = ); c = a + b*i; K = input(Type in the gain constant = ); N = input (Type in length of sequence = ); n = 1:N; x = K*exp(c*n);%Generate the sequence stem(n,real(x);%Plot the real part xlabel(Time index n);ylabel(Amplitude); title(Real part); disp(PRESS RETURN for imaginary part); pause stem(n,imag(x);%Plot the imaginary part xlabel(Time index n);ylabel(Amplitude); title(Imaginary part);% Program 2_3 % Generation of real exponential sequence % a = input(Type in argument = ); K = input(Type in the gain constant = ); N = input (Type in length of sequence = ); n = 0:N; x = K*a.n; stem(n,x); xlabel(Time index n);ylabel(Amplitude); title(alpha = ,num2str(a);% Program 2_4 % Signal Smoothing by a Moving-Average Filter % R = 50; d = rand(R,1)-0.5; m = 0:1:R-1; s = 2*m.*(0.9.m); x = s + d; plot(m,d,r-,m,s,b-,m,x,g:) xlabel(Time index n); ylabel(Amplitude) legend(dn,sn,xn); pause M = input(Number of input samples = ); b = ones(M,1)/M; y = filter(b,1,x); plot(m,s,r-,m,y,b-) legend(sn,yn); xlabel (Time index n);ylabel(Amplitude)% Program 2_5 % Illustration of Median Filtering % N = input(Median Filter Length = ); R = 50; a = rand(1,R)-0.4; b = round(a); % Generate impulse noise m = 0:R-1; s = 2*m.*(0.9.m); % Generate signal x = s + b; % Impulse noise corrupted signal y = medfilt1(x,3); % Median filtering subplot(2,1,1) stem(m,x);axis(0 50 -1 8); xlabel(n);ylabel(Amplitude); title(Impulse Noise Corrupted Signal); subplot(2,1,2) stem(m,y); xlabel(n);ylabel(Amplitude); title(Output of Median Filter);% Program 2_6 % Illustration of Convolution % a = input(Type in the first sequence = ); b = input(Type in the second sequence = ); c = conv(a, b); M = length(c)-1; n = 0:1:M; disp(output sequence =);disp(c) stem(n,c) xlabel(Time index n); ylabel(Amplitude);% Program 2_7 % Computation of Cross-correlation Sequence % x = input(Type in the reference sequence = ); y = input(Type in the second sequence = ); % Compute the correlation sequence n1 = length(y)-1; n2 = length(x)-1; r = conv(x,fliplr(y); k = (-n1):n2; stem(k,r); xlabel(Lag index); ylabel(Amplitude); v = axis; axis(-n1 n2 v(3:end);% Program 2_8 % Computation of Autocorrelation of a % Noise Corrupted Sinusoidal Sequence % N = 96; n = 1:N; x = cos(pi*0.25*n); % Generate the sinusoidal sequence d = rand(1,N) - 0.5; % Generate the noise sequence y = x + d; % Generate the noise-corrupted sinusoidal sequence r = conv(y, fliplr(y); % Compute the correlation sequence k = -28:28; stem(k, r(68:124); xlabel(Lag index); ylabel(Amplitude);% Program 3_1 % Discrete-Time Fourier Transform Computation % % Read in the desired number of frequency samples k = input(Number of frequency points = ); % Read in the numerator and denominator coefficients num = input(Numerator coefficients = ); den = input(Denominator coefficients = ); % Compute the frequency response w = 0:pi/(k-1):pi; h = freqz(num, den, w); % Plot the frequency response subplot(2,2,1) plot(w/pi,real(h);grid title(Real part) xlabel(omega/pi); ylabel(Amplitude) subplot(2,2,2) plot(w/pi,imag(h);grid title(Imaginary part) xlabel(omega/pi); ylabel(Amplitude) subplot(2,2,3) plot(w/pi,abs(h);grid title(Magnitude Spectrum) xlabel(omega/pi); ylabel(Magnitude) subplot(2,2,4) plot(w/pi,angle(h);grid title(Phase Spectrum) xlabel(omega/pi); ylabel(Phase, radians)% Program 3_2 % Generate the filter coefficients h1 = ones(1,5)/5; h2 = ones(1,14)/14; % Compute the frequency responses H1,w = freqz(h1, 1, 256); H2,w = freqz(h2, 1, 256); % Compute and plot the magnitude responses m1 = abs(H1); m2 = abs(H2); plot(w/pi,m1,r-,w/pi,m2,b-); ylabel(Magnitude); xlabel(omega/pi); legend(M=5,M=14); pause % Compute and plot the phase responses ph1 = angle(H1)*180/pi; ph2 = angle(H2)*180/pi; plot(w/pi,ph1,w/pi,ph2); ylabel(Phase, degrees);xlabel(omega/pi); legend(M=5,M=14);% Program 3_3 % Set up the filter coefficients b = -6.76195 13.456335 -6.76195; % Set initial conditions to zero values zi = 0 0; % Generate the two sinusoidal sequences n = 0:99; x1 = cos(0.1*n); x2 = cos(0.4*n); % Generate the filter output sequence y = filter(b, 1, x1+x2, zi); % Plot the input and the output sequences plot(n,y,r-,n,x2,b-,n,x1,g-.);grid axis(0 100 -1.2 4); ylabel(Amplitude); xlabel(Time index n); legend(yn,x2n,x1n)% Program 4_1 % 4-th Order Analog Butterworth Lowpass Filter Design % format long % Determine zeros and poles z,p,k = buttap(4); disp(Poles are at);disp(p); % Determine transfer function coefficients pz, pp = zp2tf(z, p, k); % Print coefficients in descending powers of s disp(Numerator polynomial); disp(pz) disp(Denominator polynomial); disp(real(pp) omega = 0: 0.01: 5; % Compute and plot frequency response h = freqs(pz,pp,omega); plot (omega,20*log10(abs(h);grid xlabel(Normalized frequency); ylabel(Gain, dB);% Program 4_2 % Program to Design Butterworth Lowpass Filter % % Type in the filter order and passband edge frequency N = input(Type in filter order = ); Wn = input(3-dB cutoff angular frequency = ); % Determine the transfer function num,den = butter(N,Wn,s); % Compute and plot the frequency response omega = 0: 200: 12000*pi; h = freqs(num,den,omega); plot (omega/(2*pi),20*log10(abs(h); xlabel(Frequency, Hz); ylabel(Gain, dB);% Program 4_3 % Program to Design Type 1 Chebyshev Lowpass Filter % % Read in the filter order, passband edge frequency % and passband ripple in dB N = input(Order = ); Fp = input(Passband edge frequency in Hz = ); Rp = input(Passband ripple in dB = ); % Determine the coefficients of the transfer function num,den = cheby1(N,Rp,2*pi*Fp,s); % Compute and plot the frequency response omega = 0: 200: 12000*pi; h = freqs(num,den,omega); plot (omega/(2*pi),20*log10(abs(h); xlabel(Frequency, Hz); ylabel(Gain, dB);% Program 4_4 % Program to Design Elliptic Lowpass Filter % % Read in the filter order, passband edge frequency, % passband ripple in dB and minimum stopband % attenuation in dB N = input(Order = ); Fp = input(Passband edge frequency in Hz = ); Rp = input(Passband ripple in dB = ); Rs = input(Minimum stopband attenuation in dB = ); % Determine the coefficients of the transfer function num,den = ellip(N,Rp,Rs,2*pi*Fp,s); % Compute and plot the frequency response omega = 0: 200: 12000*pi; h = freqs(num,den,omega); plot (omega/(2*pi),20*log10(abs(h); xlabel(Frequency, Hz); ylabel(Gain, dB);function y = circonv(x1, x2) % Develops a sequence y obtained by the circular % convolution of two equal-length sequences x1 and x2 L1 = length(x1); L2 = length(x2); if L1 = L2, error(Sequences of unequal length), end y = zeros(1,L1); x2tr = x2(1) x2(L2:-1:2); for k = 1:L1, sh = circshift(x2tr, k-1); h = x1.*sh; disp(sh); y(k) = sum(h); endfunction Y = haar_1D(X) % Computes the Haar transform of the input vector X % The length of X must be a power-of-2 % Recursively builds the Haar matrix H v = log2(length(X) - 1; H = 1 1;1 -1; for m = 1:v, A = kron(H,1 1); 2(m/2).*kron(eye(2m),1 -1); H = A; end Y = H*double(X(:);% Program 5_1 % Illustration of DFT Computation % % Read in the length N of sequence and the desired % length M of the DFT N = input(Type in the length of the sequence = ); M = input(Type in the length of the DFT = ); % Generate the length-N time-domain sequence u = ones(1,N); % Compute its M-point DFT U = fft(u,M); % Plot the time-domain sequence and its DFT t = 0:1:N-1; stem(t,u) title(Original time-domain sequence) xlabel(Time index n); ylabel(Amplitude) pause subplot(2,1,1) k = 0:1:M-1; stem(k,abs(U) title(Magnitude of the DFT samples) xlabel(Frequency index k); ylabel(Magnitude) subplot(2,1,2) stem(k,angle(U) title(Phase of the DFT samples) xlabel(Frequency index k); ylabel(Phase)% Program 5_2 % Illustration of IDFT Computation % % Read in the length K of the DFT and the desired % length N of the IDFT K = input(Type in the length of the DFT = ); N = input(Type in the length of the IDFT = ); % Generate the length-K DFT sequence k = 0:K-1; V = k/K; % Compute its N-point IDFT v = ifft(V,N); % Plot the DFT and its IDFT stem(k,V) xlabel(Frequency index k); ylabel(Amplitude) title(Original DFT samples) pause subplot(2,1,1) n = 0:N-1; stem(n,real(v) title(Real part of the time-domain samples) xlabel(Time index n); ylabel(Amplitude) subplot(2,1,2) stem(n,imag(v) title(Imaginary part of the time-domain samples) xlabel(Time index n); ylabel(Amplitude)% Program 5_3 % Numerical Computation of Fourier transform Using DFT k = 0:15; w = 0:511; x = cos(2*pi*k*3/16);% Generate the length-16 sinusoidal sequence X = fft(x); % Compute its 16-point DFT XE = fft(x,512); % Compute its 512-point DFT % Plot the frequency response and the 16-point DFT samples plot(k/16,abs(X),o, w/512,abs(XE) xlabel(omega/pi); ylabel(Magnitude)% Program 5_4 % Linear Convolution Via the DFT % % Read in the two sequences x = input(Type in the first sequence = ); h = input(Type in the second sequence = ); % Determine the length of the result of convolution L = length(x)+length(h)-1; % Compute the DFTs by zero-padding XE = fft(x,L); HE = fft(h,L); % Determine the IDFT of the product y1 = ifft(XE.*HE); % Plot the sequence generated by DFT-based convolution and % the error from direct linear convolution n = 0:L-1; subplot(2,1,1) stem(n,y1) xlabel(Time index n);ylabel(Amplitude); title(Result of DFT-based linear convolution) y2 = conv(x,h); error = y1-y2; subplot(2,1,2) stem(n,error) xlabel(Time index n);ylabel(Amplitude) title(Error sequence)% Program 5_5 % Illustration of Overlap-Add Method % % Generate the noise sequence R = 64; d = rand(R,1)-0.5; % Generate the uncorrupted sequence and add noise k = 0:R-1; s = 2*k.*(0.9.k); x = s + d; % Read in the length of the moving average filter M = input(Length of moving average filter = ); % Generate the moving average filter coefficients h = ones(1,M)/M; % Perform the overlap-add filtering operation y = fftfilt(h,x,4); % Plot the results plot(k,s,r-,k,y,b-) xlabel(Time index n);ylabel(Amplitude) legend(sn,yn)function Factors = factorize(polyn) format long; Factors = ; % Use threshold of 1e-8 instead of 0 to account for % precision effects THRESH = 1e-8; % proots = roots(polyn); % get the zeroes of the polynomial len = length(proots); % get the number of zeroes % while(len 1) if(abs(imag(proots(1) THRESH) % if the zero is a real zero fac = 1 -real(proots(1); % construct the factor with proots(1) as zero Factors = Factors;fac 0; else % if the zero has imaginary part get all zeroes whose % imag part is -ve of imaginary part of proots(1) negimag = imag(proots)+imag(proots(1); % get all zeroes which have same real part as proot(1) samereal = real(proots)-real(proots(1); %find the complex conjugate zero index = find(abs(negimag) THRESH & abs(samereal)THRESH); if(index) % if the complex conjugate exists fac = 1 -2*real(proots(1) abs(proots(1)2; %form 2nd order factor Factors = Factors;fac; else % if the complex conjugate does not exist fac = 1 -proots(1); Factors = Factors;fac 0; end end polyn = deconv(polyn,fac); %deconvolve the 1st/2nd order factor from polyn proots = roots(polyn); %determine the new zeros len = length(polyn); %determine the number of zeros end% Program 6_1 % Determination of the Factored Form % of a Rational z-Transform % num = input(Type in the numerator coefficients = ); den = input(Type in the denominator coefficients = ); K = num(1)/den(1); Numfactors = factorize(num) Denfactors = factorize(den) disp(Numerator factors);disp(Numfactors); disp(Denominator factors);disp(Denfactors); disp(Gain constant);disp(K); zplane(num,den)% Program 6_2 % Determination of the Rational z-Transform % from its Poles and Zeros % format long zr = input(Type in the zeros as a row vector = ); pr = input(Type in the poles as a row vector = ); % Transpose zero and pole row vectors z = zr; p = pr; k = input(Type in the gain constant = ); num, den = zp2tf(z, p, k); disp(Numerator polynomial coefficients); disp(num); disp(Denominator polynomial coefficients); disp(den);%Program 6_3 % Partial-Fraction Expansion of Rational z-Transform % num = input(Type in numerator coefficients = ); den = input(Type in denominator coefficients = ); r,p,k = residuez(num,den); disp(Residues);disp(r) disp(Poles);disp(p) disp(Constants);disp(k)% Program 6_4 % Partial-Fraction Expansion to Rational z-Transform % r = input(Type in the residues = ); p = input(Type in the poles = ); k = input(Type in the constants = ); num, den = residuez(r,p,k); disp(Numerator polynomial coefficients); disp(num) disp(Denominator polynomial coefficients); disp(den)% Program 6_5 % Power Series Expansion of a Rational z-Transform % % Read in the number of inverse z-transform coefficients to be computed L = input(Type in the length of output vector = ); % Read in the numerator and denominator coefficients of % the z-transform num = input(Type in the numerator coefficients = ); den = input(Type in the denominator coefficients = ); % Compute the desired number of inverse transform coefficients y,t = impz(num,den,L); disp(Coefficients of the power series expansion); disp(y)% Program 6_6 % Power Series Expansion of a Rational z-Transform % % Read in the number of inverse z-transform coefficients % to be computed N = input(Type in the length of output vector = ); % Read in the numerator and denominator coefficients of % the z-transform num = input(Type in the numerator coefficients = ); den = input(Type in the denominator coefficients = ); % Compute the desired number of inverse transform % coefficients x = 1 zeros(1, N-1); y = filter(num, den, x); disp(Coefficients of the power series expansion); disp(y)%Program_6_7 % Stability Test % num = input(Type in the numerator vector = ); den = input(Type in the denominator vector = ); N = max(length(num),length(den); x = 1; y0 = 0; S = 0;zi = zeros(1,N-1); for n = 1:1000 y,zf = filter(num,den,x,zi); if abs(y) 0.000001, break, end x = 0; S = S + abs(y); y0 = y;zi = zf; end if n 1000 disp(Stable Transfer Function); else disp(Unstable Transfer Function); end% Program 7_1 % Illustration of Deconvolution % Y = input(Type in the convolved sequence = ); H = input(Type in the convolving sequence = ); X,R = deconv(Y,H); disp(Sequence xn);disp(X); disp(Remainder Sequence rn);disp(R);% Program 7_2 % Stability Test of a Rational Transfer Function % % Read in the polynomial coefficients den = input(Type in the denominator coefficients =); % Generate the stability test parameters k = poly2rc(den); knew = fliplr(k); disp(The stability test parameters are);disp(knew); stable = all(abs(k) 1)% Program 8_2 % Factorization of a Rational IIR Transfer Function % format short num = input(Numerator coefficients = ); den = input(Denominator coefficients = ); Numfactors = factorize(num); Denfactors = factorize(den); disp(Numerator Factors),disp(Numfactors) disp(Denominator Factors),disp(Denfactors)% Program 8_3 % Parallel Realizations of an IIR Transfer Function % num = input(Numerator coefficient vector = ); den = input(Denominator coefficient vector = ); r1,p1,k1 = residuez(num,den); r2,p2,k2 = residue(num,den); disp(Parallel Form I) disp(Residues are);disp(r1); disp(Poles are at);disp(p1); disp(Constant value);disp(k1); disp(Parallel Form II) disp(Residues are);disp(r2); disp(Poles are at);disp(p2); disp(Constant value);disp(k2);% Program 8_4 % Cascaded Lattice Realization of an % Allpass Transfer Function % format long den = input(Denominator coefficient vector = ); k = poly2rc(den); knew = fliplr(k); disp(The lattice section multipliers are);disp(knew);% Program 8_5