《DSP using MATLAB》Problem 5.9

代码：

%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%%            Output Info about this m-file
fprintf('
***********************************************************
');
fprintf('        <DSP using MATLAB> Problem 5.9 

');

banner();
%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% -------------------------------------------------------
%            1  x(n) = 2*exp(-0.9*n)          n=[0:5]
%                    = 2*exp(0.9*(11-n))      n=[6:10]
%                    = 0                      else
% -------------------------------------------------------
N = 11; 
nn1 = [-5:5];
xx1 = 2*exp(-0.9*abs(nn1));

figure('NumberTitle', 'off', 'Name', 'P5.9 x(n)')
set(gcf,'Color','white'); 
%subplot(2,1,1); 
stem(nn1, xx1);
xlabel('n'); ylabel('x(n)');
title('x(n) sequence');  grid on;


n1 = [0 : 5]; n2 = [6:N-1];
x1 = 2*exp(-0.9*n1); x2 = 2*exp(-0.9*(11-n2));
x = [x1 x2 ];                              % sequence 
n = [n1  n2];
x3 = [x, zeros(1,9)];                      % padding 9 zeros
n3 = [0 : length(x3)-1];

figure('NumberTitle', 'off', 'Name', 'P5.9 x(n) and x3(n)')
set(gcf,'Color','white'); 
subplot(2,1,1); stem(n, x);
xlabel('n'); ylabel('x(n)');
title('x(n) sequence');  grid on;
subplot(2,1,2); stem(n3, x3);
xlabel('n'); ylabel('x1(n)');
title('x(n) padding zeros');  grid on;


%% =============================================================================
%%                 DTFT X(w) of xn sequence, w=[0:2pi], 
%% =============================================================================
MM = 500;
[Xw_DTFT, w] = dtft1(x, n, MM);

 magXw_DTFT = abs(Xw_DTFT); angXw_DTFT = angle(Xw_DTFT)/pi; 
realXw_DTFT = real(Xw_DTFT); imagXw_DTFT = imag(Xw_DTFT);


%% --------------------------------------------------------------
%%        START X_DTFT's  mag ang real imag
%% --------------------------------------------------------------
figure('NumberTitle', 'off', 'Name', 'P5.9 X(w) DTFT of x(n)');
set(gcf,'Color','white'); 
subplot(2,2,1); plot(w/pi,magXw_DTFT); grid on;  % axis([-2,2,0,15]); 
title('Magnitude Part');
xlabel('frequency in pi units'); ylabel('Magnitude  |X\_DTFT|'); 
subplot(2,2,3); plot(w/pi, angXw_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, realXw_DTFT); grid on;
title('Real Part');
xlabel('frequency in pi units'); ylabel('Real');
subplot('2,2,4'); plot(w/pi, imagXw_DTFT); grid on;
title('Imaginary Part');
xlabel('frequency in pi units'); ylabel('Imaginary');
%% --------------------------------------------------------------
%%        END X_DTFT's  mag ang real imag
%% --------------------------------------------------------------


%% ------------------------------------------------------------------
%%                 DFT(k) of xn sequence, k=[0:N-1]
%%                      w=2pi*k/N       k=Nw/(2pi)
%% ------------------------------------------------------------------
N1 = length(x3);
k1 = [0 : N1-1];
%k2 = [-N : N-1];
%k3 = [-N/2 : N/2];
Xk_DFT = dft(x3, N1);                                         % DFT
    magXk_DFT = abs( [ Xk_DFT ] );                          % DFT magnitude
    angXk_DFT = angle( [Xk_DFT] )/pi;                       % DFT angle
   realXk_DFT = real(Xk_DFT); imagXk_DFT = imag(Xk_DFT);

figure('NumberTitle', 'off', 'Name', 'P5.9 DFT(k) of x(n)')
set(gcf,'Color','white'); 
subplot(2,1,1); stem(k1, magXk_DFT); hold on; plot(N1*w/(2*pi), magXw_DTFT,'r--'); hold off;
%axis([-N/2, N/2, -0.5, 50.5]);
xlabel('k'); ylabel('magnitude(k)');
title('DFT magnitude of x(n), N=20');  grid on;
subplot(2,1,2); stem(k1, angXk_DFT);  hold on; plot(N1*w/(2*pi), angXw_DTFT,'r--'); hold off;
%axis([-N/2, N/2, -0.5, 50.5]);
xlabel('k'); ylabel('angle(k)');
title('DFT angle of x(n), N=20');  grid on;

　　运行结果：

       原始序列的n定义在-5到5，有负数的域。而N点序列的DFT定义在[0， N-1]中，所以需要将n进行周期N延拓，再取主值区间。

          这里在序列后补零9个，增加频率域采样点数。

        下图中，序列的DTFT（红色虚线）和DFT（蓝色棒状图）。