数字信号处理实验（一）——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