第1页 / 共35页
第2页 / 共35页
第3页 / 共35页
第4页 / 共35页
第5页 / 共35页
第6页 / 共35页
第7页 / 共35页
第8页 / 共35页
数字信号处理实践方法(第二版)例子程序(Matlab版).doc
MATLAB m-files listings – by Chapter (Digital Sign
Chapter 1 Introduction
Chapter 2 Analog I/O interface for real-time DSP s
Chapter 3 Discrete transform
Chapter 6 A framework for digital filter design
Chapter 7 Finite impulse response (FIR) filter de
Chapter 9 Multirate digital signal processing
Chapter 11 Spectrum estimation and analysis
MATLAB m-files listings – by Chapter (Digital Signal Processing – A Practical
Approach, E C Ifeachor and B W Jervis, Pearson/Prentice Hall, 2002; ISBN 81-
7808 609 3)
Chapter 1
Introduction
N/A
Chapter 2 Analog I/O interface for real-time DSP systems
N/A
%1 - forward DFT, -1 - inverse DFT
% form complex numbers
y = x*dftmtx(length(x)) ;
%compute DFT
y = x*conj(dftmtx(length(x)))/length(x);
Chapter 3 Discrete transform
function DFTD
clear all;
% Program to compute DFT coefficients directly
% (Program 3c.1, p170; program name: dftd.m)
%
direction = -1;
in = fopen('datain.dat','r');
x = fscanf(in,'%g %g',[2,inf]);
fclose(in);
x = x(1,:)+x(2,:)*i;
if direction==1
else
end
% Save/Print the results
out=fopen('dataout.dat','w');
fprintf(out,'%g %g\n',[real(y); imag(y)]);
fclose(out);
subplot(2,1,1),plot(1:length(x),x); title('Input Signal');
subplot(2,1,2),plot(1:length(y),y); title('Output Signal');
============================================================
function DFTF
%
% (Program 3c.2, p171; program name: dftf.m)
%
clear all;
direction = -1;
in = fopen('dataout.dat','r');
x = fscanf(in,'%g %g',[2,inf]);
fclose(in);
x = x(1,:)+x(2,:)*i;
if direction==1
Program to compute DFT coefficients using DIT FFT
%compute IDFT
%1 - forward DFT, -1 - inverse DFT
% form complex numbers
1
% compute FFT
% compute IFFT
y=fft(x,length(x))
y=ifft(x,length(x))
else
end
% Save/Print the results
out=fopen('dataout.dat','w');
fprintf(out,'%g %g\n',[real(y); imag(y)]);
fclose(out);
subplot(2,1,1),plot(1:length(x),x); title('Input Signal');
subplot(2,1,2),plot(1:length(y),y); title('Output Signal');
=======================================================
%
% 4-point DFT (sdft1.m)
% A simple m-file script illustrating direct 4-point DFT computation.
% input data: x(0)=1, x(1)=5, x(2)=4, x(3)=2.
%
N=4;j=sqrt(-1);
x0=1; x1=5; x2=4; x3=2;
X0=x0+x1+x2+x3;
X11=x1*exp(-j*2*pi/N);
X12=x2*exp(-j*2*pi*2/N);
X13=x3*exp(-j*2*pi*3/N);
X1a=x0+X11+X12+X13;
X1=x0+x1*exp(-j*2*pi/N)+x2*exp(-j*2*pi*2/N)+x3*exp(-j*2*pi*3/N);
X2=x0+x1*exp(-j*2*pi*2/N)+x2*exp(-j*2*pi*2*2/N)+x3*exp(-j*2*pi*2*3/N);
X3=x0+x1*exp(-j*2*pi*3/N)+x2*exp(-j*2*pi*3*2/N)+x3*exp(-j*2*pi*3*3/N);
X11
X12
X13
X1a
X0
X1
X2
X3
========================================================
%
% 4-point DFT (sdft2.m)
% A simple m-file script illustrating direct 4-point DFT computation.
% input data: x(0)=1, x(1)=0.5, x(2)=0, x(3)=0.
%
N=4;j=sqrt(-1);
x0=1; x1=0.5; x2=0; x3=0;
X0=x0+x1+x2+x3;
X11=x1*exp(-j*2*pi/N);
X12=x2*exp(-j*2*pi*2/N);
X13=x3*exp(-j*2*pi*3/N);
X1a=x0+X11+X12+X13;
X1=x0+x1*exp(-j*2*pi/N)+x2*exp(-j*2*pi*2/N)+x3*exp(-j*2*pi*3/N);
X2=x0+x1*exp(-j*2*pi*2/N)+x2*exp(-j*2*pi*2*2/N)+x3*exp(-j*2*pi*2*3/N);
X3=x0+x1*exp(-j*2*pi*3/N)+x2*exp(-j*2*pi*3*2/N)+x3*exp(-j*2*pi*3*3/N);
X11
X12
X13
2
X1a
X0
X1
X2
X3
================================================================
%
% 4-point FFT (sfft2.m)
% A simple m-file script illustrating direct 4-point FFT decimation-in-
time
% computation.
% Input data: x(0)=1, x(1)=0.5, x(2)=0, x(3)=0.
%
N=4;j=sqrt(-1);
x0=1; x1=0.5; x2=0; x3=0;
W0=1; W1=-j;
a11=x0+W0*x2;b11=x0-W0*x2;
a12=x1+W0*x3; b12=x1-W0*x3;
X0=a11+W0*a12; X2=a11-W0*a12;
X1=b11+W1*b12; X3=b11-W1*b12;
a11
b11
a12
b12
X0
X1
X2
X3
========================================================
%
%
% 8-point DFT (sfft8.m)
% A simple m-file script illustrating direct 8-point DFT computation.
% Input data: x(0)=4, x(1)=2, x(2)=1, x(3)=4, x(4)=6, x(5)=3, x(6)=5,
x(7)=2.
%
j=sqrt(-1);
x0=4; x1=2; x2=1; x3=4; x4=6; x5=3; x6=5; x7=2;
W0=1; W1=sqrt(2)/2 - j*sqrt(2)/2; W2=-j; W3=-sqrt(2)/2 - j*sqrt(2)/2;
%
A11=x0+W0*x4;B11=x0-W0*x4;
A12=x2+W0*x6;B12=x2-W0*x6;
A13=x1+W0*x5;B13=x1-W0*x5;
A14=x3+W0*x7;B14=x3-W0*x7;
% stage 2
A21=A11+W0*A12;B21=A11-W0*A12;
A22=B11+W2*B12;B22=B11-W2*B12;
A23=A13+W0*A14;B23=A13-W0*A14;
A24=B13+W2*B14;B24=B13-W2*B14;
% stage 3
X0=A21+W0*A23;X4=A21-W0*A23;
X1=A22+W1*A24;X5=A22-W1*A24;
X2=B21+W2*B23;X6=B21-W2*B23;
stage 1
3
X3=B22+W3*B24;X7=B22-W3*B24;
A11
B11
A12
B12
A13
B13
A14
B14
A21
A22
B21
B22
A23
A24
B23
B24
X0
X1
X2
X3
X4
X5
X6
X7
=================================================================
%
% 4-point inverse FFT (sifft4.m)
% A simple m-file script illustrating direct 4-point inverse DFT
computation.
% Input data: x(0)=1.5 +0j, x(1)=1-0.5j, x(2)=0.5+0j, x(3)=1+0.5j
%
N=4;j=sqrt(-1);
x0=1.5 + 0j; x1=1 - 0.5j; x2=0.5 + 0j; x3=1+0.5j;
W0=1; W1=j;
a11=x0+W0*x2;b11=x0-W0*x2;
a12=x1+W0*x3; b12=x1-W0*x3;
X0=(a11+W0*a12)/N; X2=(a11-W0*a12)/N;
X1=(b11+W1*b12)/N; X3=(b11-W1*b12)/N;
a11
b11
a12
b12
X0
X1
X2
X3
x0
x1
x2
x3
======================================================
4
% number of power series points
Chapter 4 z-transform and its applications in signal processing.
%
% m-file to compute the first 5 values of the
% inverse z-transform using the power series method
% (Program 4D.1, p235; program name: prog4d1.m)
%
n = 5;
N1 = [1 -1.122346 1]; D1 = [1 -1.433509 0.85811];
N2 = [1 1.474597 1]; D2 = [1 -1.293601 0.556929];
N3 = [1 1 0]; D3 = [1 -0.612159 0];
B = [N1; N2; N3]; A = [D1; D2; D3];
[b,a] = sos2tf([B A]);
b = [b zeros(1,n-1)];
[h,r] = deconv(b,a); %perform long division
disp(h);
============================================================
%
% m-file for finding the partial fraction expansion of a z-transform
% (Program 4D.2, p237; program name: prog4d2.m)
%
N1=[1 -1.122346 1];
N2=[1 -0.437833 1];
N3=[1 1 0];
D1=[1 -1.433509 0.85811];
D2=[1 -1.293601 0.556929];
D3=[1 -0.612159 0];
sos=[N1 D1; N2 D2; N3 D3];
[b, a] = sos2tf(sos);
[r, p, k]=residue(b, a)
================================================
%
% A script illustrating cascade to parallel realization
% (Program 4B.3, p237; program name: prog4d3.m)
%
nstage=2;
N1 =
N2 =
D1 =
D2 =
sos = [N1 D1; N2 N2];
[b, a] = sos2tf(sos);
[c, p, k] = residue(b, a);
m = length(b);
b0 = b(m)/a(m);
j=1;
for i=1:nstage
[1 0.481199 1];
[1 1.474597 1];
[1 0.052921 0.83173];
[1 -0.304609 0.238865];
bk(j)=c(j)+c(j+1);
bk(j+1)=-(c(j)*p(j+1)+c(j+1)*p(j));
ak(j)=-(p(j)+p(j+1));
ak(j+1)=p(j)*p(j+1);
j=j+2;
end
b0
5
ak
bk
c
p
k
===============================================================
%
% A simple script illustrating basics of cascade to parallel
% conversion for a 4th order filter (cprealization.m)
%
b=[1, -2, 1, 0, 0];
a=[1, -0.41421, 0.08579, 0.292895, 0.5];
[c, p, k] = residuez(b, a);
c
p
k
b0 = b(3)/a(5)
b01=c(1)+c(2)
b11=-(c(1)*p(2)+c(2)*p(1))
a11=-(p(1)+p(2))
a12=p(1)*p(2)
b02=c(3)+c(4)
b12=-(c(3)*p(4)+c(4)*p(3))
a12=-(p(3)+p(4))
a22=p(3)*p(4)
=================================================================
%
% A simple script illustrating inverse z transform by power series
codes
% File-name: pseries.m
[1 0.481199 1];
b1 =
[1 1.474597 1];
b2 =
a1 =
[1 0.052921 0.83173];
a2 =
[1 -0.304609 0.238865];
sos = [b1 a1; b2 a2];
[b, a] = sos2tf(sos);
m = length(b);
r = b;
for n=1:m
[h(n), r] = deconv(r, a);
h
r
r = [r(2:m) 0];
end
=================================================
function IZT
%program IZT (izt.m) is for:
%(1) computing the inverse z-transform via the power series or partial
%
%(2) converting a transfer function, H(z), in cascade form to an
equivalent
%
%
transfer function in parallel, via partial fraction expansion. The
basic building block is the second order biquad
fraction expansion method
6
IZT_Mode = 1;
%1 - use power series to perform the IZT
%2 - use partial fraction/residue to perform the IZT
% number of power series points
%3 - cascade to parallel conversion
n = 5;
b1 = [1 0.481199 1]; a1 = [1 0.052921 0.83173];
b2 = [1 1.474597 1]; a2 = [1 -0.304609 0.238865];
B = [b1; b2];
[b,a] = sos2tf([B A]);
if IZT_Mode==1
%compute the IZT by power series
A = [a1; a2];
b = [b zeros(1,n-1)];
[h,r] = deconv(b,a); %perform long division
disp('The result of the inverse Z-transform by power series is:');
disp(h);
else
[res,poles,rem] = residuez(b,a);
disp('The poles of the transform function are:'); disp(poles');
disp('The partial fraction coefficients are:'); disp(res');
if IZT_Mode==3
i = 1 ;
for j=1:size(B,1)
[b,a] = residuez(res(i:i+1),poles(i:i+1), 0);
fprintf('Numerator and Denominator for
stage%d:\n',j);disp(b);disp(a);
i = i + 2;
end
end
end
=========================================================
%
% m-file to illustrate the computation of frequency response
% of an IIR filter using FFT (freqrespex1.m).
%
Fs=500;
b1=[1 -1.6180 1]; b2=[];
coefficients
a1=[1 -1.5161 0.878]; a2=[];
B=[b1; b2]; A=[a1; a2];
[b,a]=sos2tf([B A]);
n=256;
dx=0;
if dx
else
% sampling frequency
% numerator/denominator filter
% number of frequency points
freqz(b,a,n,Fs);
f=(0:n-1)*(Fs/2)/n;
freqz(b,a,f,Fs);
end
=========================================================
% Frequency response evaluation by FFT
% frequency response by direct evaluation.
7
Chapter 5 Correlation and Convolution
function Correlation
%Program to compute normalised/unnormalised auto- or crosscorrelation
crosscorrelation = 1 ;
normalized = 0 ;
if crosscorrelation == 0
R11 = xcorr(x,'biased')
R11 = xcorr(x,'coeff')
%for autocorrelation
% unnormalized autocorrelation
% normalized autocorrelation
x = [1 3 2 -1 -2];
if normalized==0
else
end
x1 = [4 2 -1 3 -2 -6 -5 4 5];
x2 = [-4 1 3 7 4 -2 -8 -2 1];
if normalized==0
R12 = xcorr(x1,x2,'biased')
else
R12 = xcorr(x1,x2,'coeff') % normalized crosscorrelation
end
%for crosscorrelation
%for crosscorrelation
% unnormalized crosscorrelation
else
else
end
end
if crosscorrelation
subplot(3,1,1), plot(1:length(x1),x1), ylabel('x1(n)');
subplot(3,1,2), plot(1:length(x2),x2), ylabel('x2(n)');
subplot(3,1,3), plot(1:length(R12),R12), ylabel('R12');
subplot(2,1,1), plot(1:length(x),x), ylabel('x(n)');
subplot(2,1,2), plot(1:length(R11),R11), ylabel('R11');
Chapter 6 A framework for digital filter design
N/A
8
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
