《DSP using MATLAB》Problem 5.11

代码：

%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%%            Output Info about this m-file

fprintf('\n***********************************************************\n');

fprintf('        <DSP using MATLAB> Problem 5.11 \n\n');

banner();

%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% --------------------------------------------------------------

%              1  h(n) = [0.9^n]   n=[0:5]

% --------------------------------------------------------------

nn1 = [0:5];

hh1 = 0.9 .^ nn1;

NN1 = length(hh1);                              % length is 6

%m = mod_1(nn1, NN1);

h = [hh1 zeros(1, 2)];                        % padding zeros

n = [nn1 max(nn1)+1:max(nn1)+2];

figure('NumberTitle', 'off', 'Name', 'P5.11.1 Causal Real LSI h(n)')

set(gcf,'Color','white');

subplot(2,1,1); stem(nn1, hh1);

xlabel('n'); ylabel('h(n)');

title('h(n)=0.9^n ori sequence');  grid on;

subplot(2,1,2); stem(n, h);

xlabel('n'); ylabel('h(n)');

title('h(n) padding zeros');  grid on;

%% =============================================================================

%%                 DTFT H(w) of hn sequence, w=[0:2pi],

%% =============================================================================

MM = 500;

[Hw_DTFT, w] = dtft1(h, n, MM);

 magHw_DTFT = abs(Hw_DTFT);   angHw_DTFT = angle(Hw_DTFT)/pi;

realHw_DTFT = real(Hw_DTFT); imagHw_DTFT = imag(Hw_DTFT);

%% --------------------------------------------------------------

%%        START H_DTFT's  mag ang real imag

%% --------------------------------------------------------------

figure('NumberTitle', 'off', 'Name', 'P5.11.1 H(w) DTFT of h(n)');

set(gcf,'Color','white');

subplot(2,2,1); plot(w/pi,magHw_DTFT); grid on;  % axis([-2,2,0,15]);

title('Magnitude Part');

xlabel('frequency in \pi units'); ylabel('Magnitude  |H\_DTFT|');

subplot(2,2,3); plot(w/pi, angHw_DTFT); grid on;  % axis([-2,2,-1,1]);

title('Angle Part');

xlabel('frequency in \pi units'); ylabel('Rad \pi'); %axis([-200,200,0,2]);

subplot('2,2,2'); plot(w/pi, realHw_DTFT); grid on;

title('Real Part');

xlabel('frequency in \pi units'); ylabel('Real');

subplot('2,2,4'); plot(w/pi, imagHw_DTFT); grid on;

title('Imaginary Part');

xlabel('frequency in \pi units'); ylabel('Imaginary');

%% --------------------------------------------------------------

%%        END H_DTFT's  mag ang real imag

%% --------------------------------------------------------------

%% ------------------------------------------------------------------

%%                 DFT(k) of hn sequence, k=[0:N-1]

%%                      w=2pi*k/N       k=Nw/(2pi)

%% ------------------------------------------------------------------

N1 = length(h);

k1 = [0 : N1-1];

%k2 = [-N : N-1];

%k3 = [-N/2 : N/2];

Hk_DFT = dft(h, N1);                                        % DFT

    magHk_DFT = abs( Hk_DFT )                          % DFT magnitude

    angHk_DFT = angle( Hk_DFT )/pi                      % DFT angle

   realHk_DFT = real(Hk_DFT); imagHk_DFT = imag(Hk_DFT);

figure('NumberTitle', 'off', 'Name', 'P5.11.1 DFT(k) of h(n)')

set(gcf,'Color','white');

subplot(2,1,1); stem(k1, magHk_DFT); hold on; plot(N1*w/(2*pi), magHw_DTFT,'r--'); hold off;

%axis([-N/2, N/2, -0.5, 50.5]);

xlabel('k'); ylabel('magnitude(k)');

title('DFT magnitude of h(n), N=8');  grid on;

subplot(2,1,2); stem(k1, angHk_DFT);  hold on; plot(N1*w/(2*pi), angHw_DTFT,'r--'); hold off;

%axis([-N/2, N/2, -0.5, 50.5]);

xlabel('k'); ylabel('angle(k)');

title('DFT angle of h(n), N=8');  grid on;

　　运行结果：

只会第1小题。

8点DFT基本能反映DTFT的幅度谱，但相位谱感觉稍有欠缺。

第2小题不会做。