《DSP using MATLAB》 Problem 3.22

Posted 2020-10-21 沧海一粟

代码：

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

banner();
%% ------------------------------------------------------------------------


%% -------------------------------------------------------------------
%%           1          xa(t)=cos(20πt+θ)  through A/D
%% -------------------------------------------------------------------
   Ts = 0.05;                        % sample interval, 0.05s
   Fs = 1/Ts;                        % Fs=20Hz
%theta = 0;
%theta = pi/6;
%theta = pi/4;
%theta = pi/3;
theta = pi/2;

n1_start = 0; n1_end = 20;
      n1 = [n1_start:1:n1_end];
     nTs = n1 * Ts;              % [0, 1]s

x1 = cos(20*pi*nTs + theta * ones(1,length(n1)));            % Digital signal


M = 500;
[X1, w] = dtft1(x1, n1, M);

magX1  = abs(X1);  angX1  = angle(X1);  realX1  = real(X1);  imagX1  = imag(X1);

%% --------------------------------------------------------------------
%%              START X(w)\'s  mag ang real imag
%% --------------------------------------------------------------------
figure(\'NumberTitle\', \'off\', \'Name\', sprintf(\'Problem 3.22 X1, theta/pi = %f\', theta/pi));
set(gcf,\'Color\',\'white\'); 
subplot(2,1,1); plot(w/pi,magX1); grid on;  %axis([-1,1,0,1.05]); 
title(\'Magnitude Response\');
xlabel(\'frequency in \\pi units\'); ylabel(\'Magnitude  |H|\'); 
subplot(2,1,2); plot(w/pi, angX1/pi); grid on;  %axis([-1,1,-1.05,1.05]);
title(\'Phase Response\');
xlabel(\'frequency in \\pi units\'); ylabel(\'Radians/\\pi\');

figure(\'NumberTitle\', \'off\', \'Name\', sprintf(\'Problem 3.22 X1, theta/pi = %f\', theta/pi));
set(gcf,\'Color\',\'white\'); 
subplot(2,1,1); plot(w/pi, realX1); grid on;
title(\'Real Part\');
xlabel(\'frequency in \\pi units\'); ylabel(\'Real\');
subplot(2,1,2); plot(w/pi, imagX1); grid on;
title(\'Imaginary Part\');
xlabel(\'frequency in \\pi units\'); ylabel(\'Imaginary\');
%% -------------------------------------------------------------------
%%             END X\'s  mag ang real imag
%% -------------------------------------------------------------------

figure(\'NumberTitle\', \'off\', \'Name\', sprintf(\'Problem 3.22 xa(n), theta/pi = %f and x1(n)\', theta/pi));
na1 = 0:0.01:1;
xa1 = cos(20 * pi * na1 + theta * ones(1,length(na1)));
set(gcf, \'Color\', \'white\');
plot(1000*na1,xa1); grid on;  %axis([0,1,0,1.5]);       
title(\'x1(n) and xa(n)\');
xlabel(\'t in msec.\'); ylabel(\'xa(t)\'); hold on;
plot(1000*nTs, x1, \'o\'); hold off;


%% ------------------------------------------------------------
%%                  xa(t) reconstruction from x1(n)
%% ------------------------------------------------------------

Dt = 0.001; t = 0:Dt:1; 
xa = x1 * sinc(Fs*(ones(length(n1),1)*t - nTs\'*ones(1,length(t)))) ;

figure(\'NumberTitle\', \'off\', \'Name\', sprintf(\'Problem 3.22 Reconstructed From x1(n), theta/pi = %f\', theta/pi));
set(gcf,\'Color\',\'white\'); 
%subplot(2,1,1); 
stairs(t*1000,xa,\'r\'); grid on;  %axis([0,1,0,1.5]);       % Zero-Order-Hold
title(\'Reconstructed Signal from x1(n) using Zero-Order-Hold\');
xlabel(\'t in msec.\'); ylabel(\'xa(t)\'); hold on; 
%stem(nTs*1000, x1); gtext(\'ZOH\'); hold off;
plot(nTs*1000, x1, \'o\'); gtext(\'ZOH\'); hold off;


figure(\'NumberTitle\', \'off\', \'Name\', sprintf(\'Problem 3.22 Reconstructed From x1(n), theta/pi = %f\', theta/pi));
set(gcf,\'Color\',\'white\');
%subplot(2,1,2); 
plot(t*1000,xa,\'r\'); grid on;  %axis([0,1,0,1.5]);       % first-Order-Hold
title(\'Reconstructed Signal from x1(n) using First-Order-Hold\');
xlabel(\'t in msec.\'); ylabel(\'xa(t)\'); hold on; 
plot(nTs*1000,x1,\'o\'); gtext(\'FOH\'); hold off;



xa = spline(nTs, x1, t);
figure(\'NumberTitle\', \'off\', \'Name\', sprintf(\'Problem 3.22 Reconstructed From x1(n), theta/pi = %f\', theta/pi));
set(gcf,\'Color\',\'white\'); 
%subplot(2,1,1);
plot(1000*t, xa,\'r\'); 
xlabel(\'t in ms units\'); ylabel(\'x\');  
title(sprintf(\'Reconstructed Signal from x1(n) using Spline function\')); grid on; hold on;
plot(1000*nTs, x1,\'o\'); gtext(\'spline\');

　　运行结果：

       这里只看初相位为0的情况，原始模拟信号和采样信号（样点值圆圈标示）：

        采样信号的谱，模拟角频率20π对应的数字角频率为π，如下图所示：

        用采样信号重建原来模拟信号：

        sinc方法，stairs函数画图

        sinc方法，plot函数画图：

        cubic方法

        其他初相位的情况，这里不上图了。