数字信号处理第三章汇总

数字信号处理第三章实验程序

3.1计算离散时间傅里叶变换

% Program P3_1

% Evaluation of the DTFT

clf;

% Compute the frequency samples of the DTFT

w = -4*pi:8*pi/511:4*pi;

num = [2 1];den = [1 -0.6];

h = freqz(num, den, w);

% Plot the DTFT

subplot(2,1,1)

plot(w/pi,real(h));grid

title('Real part of H(e^{j\omega})')

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,1,2)

plot(w/pi,imag(h));grid

title('Imaginary part of H(e^{j\omega})')

xlabel('\omega /\pi');

ylabel('Amplitude');

pause

subplot(2,1,1)

plot(w/pi,abs(h));grid

title('Magnitude Spectrum |H(e^{j\omega})|')

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,1,2)

plot(w/pi,angle(h));grid

title('Phase Spectrum arg[H(e^{j\omega})]')

xlabel('\omega /\pi');

ylabel('Phase in radians');

Q3.1离散时间傅里叶变换的原始序列是H(e^jw)=(2+z^-1)/(1-0.6z^-1)。Pause的作用是暂停等待用户输入任意键后接着执行以下命令。

Q3.2

是周期函数，周期是2π。实部和幅度谱是关于y 轴对称，是偶函数；虚部和相位谱是关于原点对称，是奇函数。 Q3.3 clf; N = 512;

num = [0.7 -0.5 0.3 1]; den = [1 0.3 -0.5 0.7]; [h,w] = freqz(num, den, N); subplot(2,1,1)

plot(w/pi,real(h));grid

title('Real part of H(e^{j\omega})') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,1,2)

plot(w/pi,imag(h));grid

title('Imaginary part of H(e^{j\omega})') xlabel('\omega /\pi'); ylabel('Amplitude'); pause

subplot(2,1,1)

plot(w/pi,abs(h));grid

title('Magnitude Spectrum |H(e^{j\omega})|') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,1,2)

plot(w/pi,angle(h));grid

title('Phase Spectrum arg[H(e^{j\omega})]') xlabel('\omega /\pi');

ylabel('Phase in radians');

还是周期函数，周期是2π。相位谱的跳变的原因是：在利用反正切函数计算角度的时候，其中的一个分支出现了衰减，造成了跳变。

clf;

N = 512;

num = [0.7 -0.5 0.3 1];

den = [1 0.3 -0.5 0.7];

[h,w] = freqz(num, den, N);

subplot(2,1,1)

plot(w/pi,unwrap(angle(h)));grid

title('Phase Spectrum arg[H(e^{j\omega})]')

xlabel('\omega /\pi');

ylabel('Phase in radians');

Q3.4 修改后的程序为

clf;

w = -4*pi:8*pi/511:4*pi;

num = [1 3 5 7 9 11 13 15 17];

den = 1;

h = freqz(num, den, w);

% Plot the DTFT

subplot(2,1,1)

plot(w/pi,real(h));grid

title('Real part of H(e^{j\omega})')

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,1,2)

plot(w/pi,imag(h));grid

title('Imaginary part of H(e^{j\omega})')

xlabel('\omega /\pi');

ylabel('Amplitude');

pause

subplot(2,1,1)

plot(w/pi,abs(h));grid

title('Magnitude Spectrum |H(e^{j\omega})|')

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,1,2)

plot(w/pi,angle(h));grid

title('Phase Spectrum arg[H(e^{j\omega})]')

xlabel('\omega /\pi');

ylabel('Phase in radians');

是周期函数，周期是2π。实部和幅度谱是关于y轴对称，是偶函数；虚部和相位谱是关于原点对称，是奇函数。

Q3.5若要改为以度为单位，则将程序中的第二个图的程序改为

subplot(2,1,2)

plot(w/pi,180*angle(h)/pi);grid

title('Phase Spectrum arg[H(e^{j\omega})]')

xlabel('\omega /\pi');

ylabel('Phase in degrees');

就可以了。

3.2离散时间傅里叶变换的性质

1.时移特性

clf;

w = -pi:2*pi/255:pi;

D = 10;

num = [1 2 3 4 5 6 7 8 9];

h1 = freqz(num, 1, w);

h2 = freqz([zeros(1,D) num], 1, w);

subplot(2,2,1)

plot(w/pi,abs(h1));grid

title('Magnitude Spectrum of Original Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,2,2)

plot(w/pi,abs(h2));grid

title('Magnitude Spectrum of Time-Shifted Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,2,3)

plot(w/pi,angle(h1));grid

title('Phase Spectrum of Original Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Phase in radians');

subplot(2,2,4)

plot(w/pi,angle(h2));grid

title('Phase Spectrum of Time-Shifted Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Phase in radians');

Q3.6参数D控制时移量。

Q3.7

D=10 D=50

时移特性：信号在时域移动某个距离，则所得信号的幅度谱和原信号相同，而相位谱是原信号的相位谱再附加一个线性相移，由时移特性可以看到，信号的相位谱可以反映信号在时域中的位置信息，不同位置上的同一信号，它们具有不同的相频特性，而幅频特性相同。

Q3.8如上图所示

Q3.9改变序列长度

num = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 2425 26 27 28 29];所得的图像为

D=10 D=50

从上图中可以看出，增加序列的长度，使得幅度谱更加窄，而相位谱则更加密集和陡峭。

2.平移特性

Q3.10

clf;

w = -pi:2*pi/255:pi;

wo = 0.4*pi;

num1 = [1 3 5 7 9 11 13 15 17];

L = length(num1);

h1 = freqz(num1, 1, w);

n = 0:L-1;

num2 = exp(wo*i*n).*num1;

h2 = freqz(num2, 1, w);

subplot(2,2,1)

plot(w/pi,abs(h1));grid

title('Magnitude Spectrum of Original Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,2,2)

plot(w/pi,abs(h2));grid

title('Magnitude Spectrum of Frequency-Shifted Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,2,3)

plot(w/pi,angle(h1));grid

title('Phase Spectrum of Original Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Phase in radians');

subplot(2,2,4)

plot(w/pi,angle(h2));grid

title('Phase Spectrum of Frequency-Shifted Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Phase in radians');

Wo控制平移量。

-1

-0.500.51

0100200

300序列num1的幅度谱

ω/π

振幅

-1

-0.5

00.510100200

300频移(wo=0.8*pi)后序列的幅度谱ω/π

振幅

-1

-0.5

00.5

1

-4-202

4序列num1的相位谱

ω/π

振幅

-1

-0.5

00.5

1

-4-202

4频移(wo=0.8*pi)后序列的相位谱ω/π

振幅

-1

-0.5

00.51

010*******

400序列num2的幅度谱

ω/π

振幅

-1

-0.5

0.5

1

010*******

400频移(wo=-0.2*pi)后序列的幅度谱ω/π

振幅

-1

-0.5

00.5

1

-4-202

4序列num2的相位谱ω/π

振幅

-1

-0.5

0.5

1

-4

-2

02

4频移(wo=-0.2*pi)后序列的相位谱ω/π

振幅

Q3.11由结果图Q3.11可得出在参数wo 的控制下，离散时间傅里叶变换的幅度谱和相位谱都随着控制参数右移k 个单位（wo=k*pi ）。

k=0.4 k=-0.4

Q3.12将k 改为-0.4得到的运行结果如上图。 Q3.13改变序列长度

序列：num1=[1 3 5 7 9 11 13 15 17 19 21 23 25 27 29]

序列：num2=[11 13 15 17 19 21 23 25 27 29 31 33 35 37 39];

3.卷积性质 Q3.14 clf;

w = -pi:2*pi/255:pi; % freqency vector for evaluating DTFT x1 = [1 3 5 7 9 11 13 15 17]; x2 = [1 -2 3 -2 1]; y = conv(x1,x2);

h1 = freqz(x1, 1, w); h2 = freqz(x2, 1, w); hp = h1.*h2;

h3 = freqz(y,1,w); subplot(2,2,1)

plot(w/pi,abs(hp));grid

title('Product of Magnitude Spectra','FontSize',8)

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,2,2)

plot(w/pi,abs(h3));grid

title('Magnitude Spectrum of Convolved Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,2,3)

plot(w/pi,angle(hp));grid

title('Sum of Phase Spectra','FontSize',8)

xlabel('\omega /\pi');

ylabel('Phase in radians');

subplot(2,2,4)

plot(w/pi,angle(h3));grid

title('Phase Spectrum of Convolved Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Phase in radians');

Q3.15

分析结果图可以得出幅度谱的乘积和卷积后的幅度谱相同，相位谱的乘积和卷积后的相位谱相同。

Q3.16 x1 = [1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33];

x2 = [1 -2 3 -2 1 -5 2 -3 1];运行结果如上边第二个图所示。

4.调制性质

Q3.17

clf;

w = -pi:2*pi/255:pi;

x1 = [1 3 5 7 9 11 13 15 17];

x2 = [1 -1 1 -1 1 -1 1 -1 1];

y = x1.*x2;

h1 = freqz(x1, 1, w);

h2 = freqz(x2, 1, w);

h3 = freqz(y,1,w);

subplot(3,1,1)

plot(w/pi,abs(h1));grid

title('Magnitude Spectrum of First Sequence')

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(3,1,2)

plot(w/pi,abs(h2));grid

title('Magnitude Spectrum of Second Sequence')

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(3,1,3)

plot(w/pi,abs(h3));grid

title('Magnitude Spectrum of Product Sequence')

xlabel('\omega /\pi');

ylabel('Amplitude');

Q3.18分析图得出乘积序列的幅度谱近似等于两序列的幅度谱的和.

Q3.19将序列改变为x1 = [1 3 5 7 9 11 13 15 17 19 21 23 25 27]; x2 = [1 -1 1 -1 1 -1 1 -1 1 0 2 -4 7 -1]得到的运行结果为上右图。乘积序列的幅度谱依然近似等于两序列的幅度谱的和.

5.时间反转性质

Q3.20

clf;

w = -pi:2*pi/255:pi;

num = [1 2 3 4];

L = length(num)-1;

h1 = freqz(num, 1, w);

h2 = freqz(fliplr(num), 1, w);

h3 = exp(w*L*i).*h2;

subplot(2,2,1)

plot(w/pi,abs(h1));grid

title('Magnitude Spectrum of Original Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,2,2)

plot(w/pi,abs(h3));grid

title('Magnitude Spectrum of Time-Reversed Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Amplitude');

subplot(2,2,3)

plot(w/pi,angle(h1));grid

title('Phase Spectrum of Original Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Phase in radians');

subplot(2,2,4)

plot(w/pi,angle(h3));grid

title('Phase Spectrum of Time-Reversed Sequence','FontSize',8)

xlabel('\omega /\pi');

ylabel('Phase in radians');

Q3.21分析图得出序列的幅度谱随时间反转不发生变化，序列相位谱随时间反转而反转180。

Q3.22改变序列长度num = [1 -2 3 -4 5 -6 7 -8 ];得到的运行结果为上右，结果依然是序列的幅度谱随时间反转不发生变化，序列相位谱随时间反转而反转180。

3.5离散傅里叶变换和离散傅里叶逆变换的运算

Q3.23

clf;

N=200;

L=256;

nn = [0:N-1];

kk = [0:L-1];

xR = [0.1*(1:100) zeros(1,N-100)]; xI = [zeros(1,N)];

x = xR + i*xI;

XF = fft(x,L);

subplot(3,2,1);grid;

plot(nn,xR);grid;

title('Re\{x[n]\}');

xlabel('Time index n');

ylabel('Amplitude');

subplot(3,2,2);

plot(nn,xI);grid;

title('Im\{x[n]\}');

xlabel('Time index n');

ylabel('Amplitude');

subplot(3,2,3);

plot(kk,real(XF));grid;

title('Re\{X[k]\}');

xlabel('Frequency index k');

ylabel('Amplitude');

subplot(3,2,4);

plot(kk,imag(XF));grid;

title('Im\{X[k]\}');

xlabel('Frequency index k');

ylabel('Amplitude');

xx = ifft(XF,L);

subplot(3,2,5);

plot(kk,real(xx));grid;

title('Real part of IDFT\{X[k]\}'); xlabel('Time index n');

ylabel('Amplitude');

subplot(3,2,6);

plot(kk,imag(xx));grid;

title('Imag part of IDFT\{X[k]\}'); xlabel('Time index n');

ylabel('Amplitude');

Q3.24

clf;

N=256;

nn = [0:N-1];

ntime = [-N/2:N/2-1];

g = (0.75).^abs(ntime);

h = (-0.9).^ntime;

GF = fft(g);

HF = fft(h);

x = g + i*h;

XF = fft(x);

XFstar = conj(XF);

XFstarmod = [XFstar(1) fliplr(XFstar(2:N))]; GF2 = 0.5*(XF + XFstarmod);

HF2 = -i*0.5*(XF - XFstarmod);

abs(max(GF-GF2))

abs(max(HF-HF2))

figure(1);clf;

subplot(2,2,1);grid;

plot(nn,real(GF));grid;

title('Two N-point DFT''s');

xlabel('Frequency index k');

ylabel('Re\{G[k]\}');

subplot(2,2,2);

plot(nn,imag(GF));grid;

title('Two N-point DFT''s');

xlabel('Frequency index k');

ylabel('Im\{G[k]\}');

subplot(2,2,3);grid;

plot(nn,real(GF2));grid;

title('Single N-point DFT');

xlabel('Frequency index k');

ylabel('Re\{G[k]\}');

subplot(2,2,4);

plot(nn,imag(GF2));grid;

title('Single N-point DFT'); xlabel('Frequency index k'); ylabel('Im\{G[k]\}');

figure(2);clf;

subplot(2,2,1);grid;

plot(nn,real(HF));grid;

title('Two N-point DFT''s'); xlabel('Freq index k');

ylabel('Re\{H[k]\}');

subplot(2,2,2);

plot(nn,imag(HF));grid;

title('Two N-point DFT''s'); xlabel('Freq index k');

ylabel('Im\{H[k]\}');

subplot(2,2,3);grid;

plot(nn,real(HF2));grid;

title('Single N-point DFT'); xlabel('Freq index k');

ylabel('Re\{H[k]\}');

subplot(2,2,4);

plot(nn,imag(HF2));grid;

title('Single N-point DFT'); xlabel('Freq index k');

ylabel('Im\{H[k]\}');

Q3.25

clf;

N = 128;

TwoN = 2*N;W2N = exp(-i*pi/N); k = [0:TwoN-1];

v = (-0.7.^k);

g = downsample(v,2);

h = downsample(v,2,1);

x = g + i*h;

XF = fft(x);

XFstar = conj(XF);

XFstarmod = [XFstar(1) fliplr(XFstar(2:N))];

GF = 0.5*(XF + XFstarmod);

HF = -i*0.5*(XF - XFstarmod);

VF = [GF GF] + (W2N.^k).*[HF HF];

VF2 = fft(v);

abs(max(VF-VF2))

subplot(2,2,1);

plot(k,real(VF));grid;

title('Complex N-point DFT');

xlabel('Frequency index k');

ylabel('Re\{V[k]\}');

subplot(2,2,2);

plot(k,imag(VF));grid;

title('Complex N-point DFT');

xlabel('Frequency index k');

ylabel('Im\{V[k]\}');

subplot(2,2,3);

plot(k,real(VF2));grid;

title('Real 2N-point DFT');

xlabel('Frequency index k');

ylabel('Re\{V[k]\}');

subplot(2,2,4);

plot(k,imag(VF2));grid;

title('Real 2N-point DFT');

xlabel('Frequency index k');

ylabel('Im\{V[k]\}');

3.4离散傅里叶函数的性质

Q3.26rem(x,y),x是除y以后剩余部分。

Q3.27输入序列x循环移位留下的位置。如果M > 0,那么circshift删除左边的元素向量x

和附加他们右侧获得剩下的元素循环转移序列。如果如果M < 0,然后circshift第一次补充的x的长度,即。,最右边的长度(x)- m样品从x和附加右边的样品得到循环转移序列。Q3.28这是二元关系不等于操作符。~ = B返回值1如果A和B是不平等的值0如果A和B 都是平等的。

Q3.29输入是平等的两个向量x1和x2长度l .理解circonv是如何工作的,它是有用的定期x2的延伸。让x2p x2的无限长的周期延长。从概念上讲,常规时间逆转x2p和集x2tr 1到L等于元素的时间逆转x2p版本。元素1通过y L的输出向量然后通过x1和长度之间的内积向量sh循环变化对时间逆转向量x2tr。对于输出样例y[n],1≤n≤L、正确的循环移位是n - 1的位置。

Q3.30

clf;

M = 6;

a = [0 1 2 3 4 5 6 7 8 9];

b = circshift(a,M);

L = length(a)-1;

n = 0:L;

subplot(2,1,1);

stem(n,a);axis([0,L,min(a),max(a)]);

title('Original Sequence');

xlabel('time index n');

ylabel('a[n]');

subplot(2,1,2);

stem(n,b);axis([0,L,min(a),max(a)]);

title(['Sequence Obtained by Circularly Shifting by ',num2str(M),'Samples']); xlabel('time index n');

ylabel('b[n]')

M值决定时移量。

Q3.31

Q3.32

clf;

x = [0 2 4 6 8 10 12 14 16];

N = length(x)-1; n = 0:N;

y = circshift(x,5);

XF = fft(x);

YF = fft(y);

subplot(2,2,1);

stem(n,abs(XF));grid;

title('Magnitude of DFT of Original Sequence');

xlabel('Frequency index k');

ylabel('|X[k]|');

subplot(2,2,2);

stem(n,abs(YF));grid;

title('Magnitude of DFT of Circularly Shifted Sequence'); xlabel('Frequency index k');

ylabel('|Y[k]|');

subplot(2,2,3);

stem(n,angle(XF));grid;

title('Phase of DFT of Original Sequence');

xlabel('Frequency index k');

ylabel('arg(X[k])');

subplot(2,2,4);

stem(n,angle(YF));grid;

title('Phase of DFT of Circularly Shifted Sequence'); xlabel('Frequency index k');

ylabel('arg(Y[k])');

时移量是8.

Q3.33

Q3.34 M=5 运行结果如上右图所示。

Q3.35Length = 13 Length = 20

Q3.36

g1 = [1 2 3 4 5 6]; g2 = [1 -2 3 3 -2 1];

ycir = circonv(g1,g2);

disp('Result of circular convolution = ');disp(ycir)

G1 = fft(g1); G2 = fft(g2);

yc = real(ifft(G1.*G2));

disp('Result of IDFT of the DFT products = ');disp(yc)

运行结果

Q3.37结果如下：

Q3.38

g1 = [1 2 3 4 5];g2 = [2 2 0 1 1];

g1e = [g1 zeros(1,length(g2)-1)];

g2e = [g2 zeros(1,length(g1)-1)];

ylin = circonv(g1e,g2e);

disp('Linear convolution via circular convolution = ');disp(ylin);

y = conv(g1, g2);

disp('Direct linear convolution = ');disp(y)

结果如下：

Q3.39 g1 = [3 1 4 1 5 9 2];g2 = [1 1 1 0 0]; g1 = [5 4 3 2 1 0];g2 = [-2 1 2 3 4];

Q3.40

g1 = [1 2 3 4 5];

g2 = [2 2 0 1 1];

g1e = [g1 zeros(1,length(g2)-1)];

g2e = [g2 zeros(1,length(g1)-1)];

G1EF = fft(g1e);

G2EF = fft(g2e);

ylin = real(ifft(G1EF.*G2EF));

disp('直线线性卷积 = ' );disp(ylin)；

Q3.41

x = [1 2 4 2 6 32 6 4 2 zeros(1,247)];

x1 = [x(1) x(256:-1:2)];

xe = 0.5 *(x + x1);

XF = fft(x);

XEF = fft(xe);

clf;

k = 0:255;

subplot(2,2,1);

plot(k/128,real(XF)); grid;

ylabel('Amplitude');

title('Re(DFT\{x[n]\})');

subplot(2,2,2);

plot(k/128,imag(XF)); grid;

ylabel('Amplitude');

title('Im(DFT\{x[n]\})');

subplot(2,2,3);

plot(k/128,real(XEF)); grid;

xlabel('Time index n');ylabel('Amplitude');

title('Re(DFT\{x_{e}[n]\})');

subplot(2,2,4);

plot(k/128,imag(XEF)); grid;

xlabel('Time index n');ylabel('Amplitude');

title('Im(DFT\{x_{e}[n]\})');

X1[n]是X[n]的逆转序列

Q3.42 XEF等于零的虚部在浮点精度。这个结果可以解释如下:真正的转换的一部分x[n]的转换定期甚至x[n]的一部分。因此,定期的DFT甚至x[n]的一部分，真正的一部分,正是真正的X[k]的一部分,一个虚部为零。

Q3.43仿真图为上右图。

Q3.44a和b的值相等。

x = [(1:128) (128:-1:1)];

XF = fft(x);

a = sum(x.*x)

b = round(sum(abs(XF).^2)/256)

Q3.45

x = [(1:128) (128:-1:1)];

XF = fft(x);

a = sum(x.*x)

b = round(sum(XF.*conj(XF))/256)

Q3.46

Q3.47运行结果为

Zeros:

-1.0000 + 1.4142i

-1.0000 - 1.4142i

-0.2500 + 0.6614i

-0.2500 - 0.6614i

Poles:

-8.9576

-0.2718

0.1147 + 0.2627i

0.1147 - 0.2627i

sos =

1.0000

2.0000

3.0000 1.0000 9.2293 2.4344 1.0000 0.5000 0.5000 1.0000 -0.2293 0.0822 k =

0.4000

Q3.48 结果为：

Q3.49