數字信號處理實驗（一）——DTFT

一、離散序列傅里葉變化——DTFT

1、DTFT公式

2、Matlab算法實現

function[X]=dtft(x,n,w,flag)

%計算離散時間付里葉變換
%[X]=dtft(x,n,w)
%X=在w頻率點上的DTFT數組
%x=n點有限長度序列
%n=樣本位置向量
%w=頻率點位置向量

X = x * (exp(-j).^(n' * w));

3、DTFT一些畫圖代碼

function [] = signal_write(X,w,flag)
% X:數據
% w:頻率向量

 magX=abs(X);angX=angle(X);
 realX=real(X);imagX=imag(X);

if(flag == 1)
    figure();
    magX=abs(X);angX=angle(X);
    realX=real(X);imagX=imag(X);
    subplot(2,2,1);plot(w/pi,magX);grid
    xlabel('以pi為單位的頻率');title('幅度部分');ylabel('幅度')
    subplot(2,2,3);plot(w/pi,angX);grid
    xlabel('以pi為單位的頻率');title('相角部分');ylabel('弧度')
    subplot(2,2,2);plot(w/pi,realX);grid
    xlabel('以pi為單位的頻率');title('實部');ylabel('實部')
    subplot(2,2,4);plot(w/pi,imagX);grid
    xlabel('以pi為單位的頻率');title('虛部');ylabel('虛部')
end

if(flag == 2)
    plot(w/pi,magX);grid
    xlabel('以pi為單位的頻率');title('幅度部分');ylabel('幅度')
end


if(flag == 3)
    plot(w/pi,angX);grid
    xlabel('以pi為單位的頻率');title('相角部分');ylabel('弧度')
end

if(flag == 4)
    plot(w/pi,realX);grid
    xlabel('以pi為單位的頻率');title('實部');ylabel('實部')
end

if(flag == 5)
    plot(w/pi,imagX);grid
    xlabel('以pi為單位的頻率');title('虛部');ylabel('虛部')
end

 

 

二、求LTI系統的頻率響應H

%example2.4
clear all;close all;

b=[1];
a=[1 -0.8];

m=0;length(b)-1;
l=0:length(a)-1;

%頻率分點
K=500;
k=-2*K:1:2*K;
w=pi*k/K;

%構建分子和分母的傅里葉變換
num=b*exp(-j*m'*w); %分母
den=a*exp(-j*l'*w); %分子
h=num./den;

magH=abs(h);
angH=angle(h);

figure(1)
subplot(2,1,1);plot(w/pi,magH),grid,title('幅度部分')
subplot(2,1,2);plot(w/pi,angH),grid,title('相角部分')
n=0:100;
x=cos(0.05*pi*n);
y=filter(b,a,x);
figure(2)
subplot(2,1,1);plot(n,x),grid,title('輸入信號')
subplot(2,1,2);plot(n,y),grid,title('輸出信號')

 

三、采樣與重構

Matlab代碼

function [  ] = caiyang(Fs,N,jt,flag)
%UNTITLED3 此處顯示有關此函數的摘要
%   此處顯示詳細說明
%   Dt  模擬時間間隔：在特定精度下信號為模擬的 
%   t   模擬時刻序列 
%   n   離散時間索引
%   Ts  采樣周期 
%   Fs  采樣頻率
%   xa  在特定精度下，由離散信號逼近模擬信號
%   jt  是否需要重構
%   flag  5  第五題 
%         6  第六題

Dt=0.00005;             % 模擬時間間隔：在特定精度下信號為模擬的
Ts=1/Fs;                % 采樣周期 
n=-N:1:N;               % 離散時間索引
nTs=n*Ts;               % 序列時刻索引
t=-N*Ts:Dt:N*Ts;        % 模擬時刻序列 

%% 只是對應相應的作業、、
if flag == 5
   xa=exp(-1000*abs(t));       % 在特定精度下，由離散信號逼近模擬信號

   x1=exp(-1000*abs(nTs));     % Fs=5000 樣本/s:x1為采樣后的離散時間序列 
end

if flag == 6
   xa=sin(1000*pi*t);       % 在特定精度下，由離散信號逼近模擬信號

   x1=sin(1000*pi*nTs);     % Fs=5000 樣本/s:x1為采樣后的離散時間序列 
   
end

if flag == 7
    xa = sin(20*pi*t + pi/4);
    x1 = sin(20*pi*nTs + pi/4);
end

%%

K  = 500; % 對pi進行K等分:相當於對單位園2pi進行1000等分 
dk = pi/K; % pi 的等分步進索引 
w  = 0 : dk : pi; % 角度步進索引 
X  = x1 * exp(-j* n'*w);% 對x1序列做離散傅立葉變換 
Xr = real(X);
w  = [-fliplr(w),w(2:end)]; % 定義w負半軸 
Xr = [fliplr(Xr),Xr(2:end)]; % 由於實部偶對稱，得到Xr的負半軸 

%% 決定是否重構
if jt == 1
   figure();
   % 繪出xa
   subplot(3,1,1)
   plot(t*1000,xa);hold on         
   % 繪出x(jw)
   stem(n*Ts*1000,x1,'r.'),hold off,title('時域波形') 

   % 繪出以pi歸一化的數字頻率對應的頻域實部波形
   subplot(3,1,2);plot(w/pi,Xr);title('頻域波形')     
   
   subplot(3,1,3)
   chonggou(x1,Fs,N);
end

if jt == 0
   figure();
   % 繪出xa
   subplot(2,1,1);
   plot(t*1000,xa);hold on         
   % 繪出x(jw)
   stem(n*Ts*1000,x1,'r.'),hold off,title('時域波形') 

   % 繪出以pi歸一化的數字頻率對應的頻域實部波形
   subplot(2,1,2);plot(w/pi,Xr);title('頻域波形')      
end
   
if jt == 2
   % 繪出以pi歸一化的數字頻率對應的頻域實部波形
   plot(w/pi,Xr);title('頻域波形')    

end

if jt == 3 
   figure();
   subplot(2,1,1);
   % 繪出xa
   plot(t*1000,xa);hold on         
   % 繪出x(jw)
   plot(n*Ts*1000,x1,'r.'),hold off,title('時域波形') 
 
   xa  = x1 * sinc(Fs*(ones(length(nTs),1) * t-nTs'*ones(1,length(t))));      % 內插重構
   
   subplot(2,1,2);
   plot(t*1000,xa, 'k' ),hold on
   plot(n*Ts*1000,x1,'r.'),hold off ,title('重構波形' )
   axis([-N/Fs*1000,N/Fs*1000,min(x1),max(x1)]);
   




end

 

重構代碼：

function [  ] = chonggou(x1,Fs,N)
%UNTITLED4 此處顯示有關此函數的摘要
%   此處顯示詳細說明
%   x1 抽樣序列
%   Fs 采樣率
%   t  時間軸
%   Dt 離散間隔，模擬信號


Dt  = 0.00005;             % 模擬時間間隔：在特定精度下信號為模擬的
n   = -N:N;
nTs = n/Fs;
t   = -N/Fs:Dt:N/Fs;        % 模擬時刻序列 
xa  = x1 * sinc(Fs*(ones(length(nTs),1) * t-nTs'*ones(1,length(t))));      % 內插重構

plot(t*1000,xa, 'k' ),hold on
stem(nTs*1000,x1, 'r.' ),hold off ,title('重構波形' )
axis([-N/Fs*1000,N/Fs*1000,min(x1),max(x1)]);

end