代码:
%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ %% Output Info about this m-file fprintf('\n***********************************************************\n'); fprintf(' <DSP using MATLAB> Problem 5.38 \n\n'); banner(); %% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % ------------------------------------------------------------------------ % 1 x(n) = [2sin(0.04pi*n)+5cos(0.08pi*n)] % N=40 N=50 N=60 % % ------------------------------------------------------------------------ N = 40; nn1 = [0:N-1]; xx1 = 2*sin(0.04*pi*nn1)+5*cos(0.08*pi*nn1); NN1 = length(xx1); % length %m = mod_1(nn1, NN1); %x = [xx1 zeros(1, 0)]; % padding zeros %n = [nn1 max(nn1)+1:max(nn1)+6]; x = xx1; n = nn1; str = fprintf('P5.38.1 x(n) N=%d', N); figure('NumberTitle', 'off', 'Name', str); set(gcf,'Color','white'); subplot(2,1,1); stem(nn1, xx1); xlabel('n'); ylabel('x(n)'); title('x(n) ori sequence'); grid on; subplot(2,1,2); stem(n, x); xlabel('n'); ylabel('x(n)'); title('x(n) padding 0 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.38.1 X(w) DTFT of x(n) N=40'); 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(x); k1 = [0 : N1-1]; %k2 = [-N : N-1]; %k3 = [-N/2 : N/2]; Xk_DFT = dft(x, 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.38.1 DFT(k) of x(n) N=40') 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)'); 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)'); grid on; %% ============================================================================= %% DTFT X(w) of xn sequence, w=[-2pi:2pi], %% ============================================================================= MM = 500; [Xw_DTFT, w] = dtft2(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.38.1 X(w) DTFT of x(n) N=40'); 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=[-N:N-1] %% w=2pi*k/N k=Nw/(2pi) %% ------------------------------------------------------------------ N2 = length(x); %k1 = [0 : length(x)-1]; k2 = [-N2 : N2-1]; %k3 = [-N/2 : N/2]; Xk_DFT = dft(x, N2); % DFT magXk_DFT = abs( [ Xk_DFT Xk_DFT] ); % DFT magnitude angXk_DFT = angle( [Xk_DFT Xk_DFT] )/pi; % DFT angle realXk_DFT = real([Xk_DFT Xk_DFT]); imagXk_DFT = imag([Xk_DFT Xk_DFT]); figure('NumberTitle', 'off', 'Name', 'P5.38.1 DFT(k) of x(n) N=40') set(gcf,'Color','white'); subplot(2,1,1); stem(k2, magXk_DFT); hold on; plot(N2*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)'); grid on; subplot(2,1,2); stem(k2, angXk_DFT); hold on; plot(N2*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)'); grid on; %% ============================================================================= %% DTFT X(w) of xn sequence, w=[-pi:pi], %% ============================================================================= MM = 500; [Xw_DTFT, w] = dtft3(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.38.1 X(w) DTFT of x(n) N=40'); 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=[-N/2:N/2] %% w=2pi*k/N k=Nw/(2pi) %% ------------------------------------------------------------------ N3 = length(x); %k1 = [0 : length(x)-1]; %k2 = [-N : N-1]; k3 = [-N3/2 : N3/2]; Xk_DFT = dft(x, N3); % DFT magXk_DFT = abs( [ Xk_DFT(N3/2+1 : N3) Xk_DFT(1 : N3/2+1)] ); % DFT magnitude angXk_DFT = angle( [Xk_DFT(N3/2+1 : N3) Xk_DFT(1 : N3/2+1)] )/pi; % DFT angle % realXk_DFT = real([Xk_DFT(N/2+1 : N) Xk_DFT(1 : N/2+1)]); % imagXk_DFT = imag([Xk_DFT(N/2+1 : N) Xk_DFT(1 : N/2+1)]); figure('NumberTitle', 'off', 'Name', 'P5.38.1 DFT(k) of x(n) N=40') set(gcf,'Color','white'); subplot(2,1,1); stem(k3, magXk_DFT); hold on; plot(N3*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)'); grid on; subplot(2,1,2); stem(k3, angXk_DFT); hold on; plot(N3*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)'); grid on;
运行结果:
下面是N=50的情况,附图