1.dfdt
function X =dtft(x,n,w) %计算离散时间付里叶变换 %[X]=dtft(x,n,w) %X=在w频率点上的DTFT数组 %x=n点有限长度序列 %n=样本位置向量 %w=频率点位置向量 X=x*(exp(-j).^(n'*w));
2.idfdt
function[x]=idtft(X,n,w) %计算离散时间付里叶变换 %[X]=dtft(x,n,w) %X=在w频率点上的DTFT数组 %x=n点有限长度序列 %n=样本位置向量 %w=频率点位置向量 x=1/(2*pi)*X*(exp(j).^(w'*n));
3.sigfold
function [y, n] = sigfold(x, n) y = fliplr(x); n = -fliplr(n);
4.sigshift
function [y,n] = sigshift(x,m,n0) % 实现 y(n) = x(n-n0) % ------------------------- % [y,n] = sigshift(x,m,n0) % n = m+n0; y = x;
5.绘制dtft的幅度、相角、实部、虚部
%% (a) n=0:20; x=2*(0.8).^n; w=0:pi/500:pi; X=dtft(x,n,w) magX=abs(X); angX=angle(X); realX=real(X); imagX=imag(X); figure(1) 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('虚部')
6.dtft的性质验证
%% 3 n=0:10; x1=rand(1,length(n)); x2=rand(1,length(n)); x=x1; w=0:pi/500:pi; X1=dtft(x1,n,w); X2=dtft(x2,n,w); X=dtft(x,n,w);
%% 线性 a=0.3; b=0.7; F1=dtft(a*x1+b*x2,n,w); F2=a*X1+b*X2; magF1=abs(F1); angF1=angle(F1); magF2=abs(F2); angF2=angle(F2); figure(1) subplot(2,2,1); plot(w/pi,magF1);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,3);plot(w/pi,angF1);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度') subplot(2,2,2); plot(w/pi,magF2);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,4);plot(w/pi,angF2);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度')
%% 时移 k=5; %F1=dtft(x,n-k,w); F2=exp(-j*w*k).*X1; y2=idtft(F2,n,w); %[y1,n1]=sigshift(x,n,k); %n2=n; subplot(2,1,1);stem(n,x,'r.'),title('原始时域波形') % 绘出x(jw) subplot(2,1,2);stem(n2,y2,'r.'),title('平移后时域波形') % 绘出x(jw)
%% 频移 w0=pi/3; F1=dtft(x.*exp(j*w0*n),n,w); F2=dtft(x,n,w); magF1=abs(F1); angF1=angle(F1); magF2=abs(F2); angF2=angle(F2); figure(3) subplot(2,2,1); plot(w/pi,magF1);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,3);plot(w/pi,angF1);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度') subplot(2,2,2); plot(w/pi,magF2);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,4);plot(w/pi,angF2);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度')
%% 共轭 F1=X1; F2=conj(dtft(x,n,-w)); magF1=abs(F1); angF1=angle(F1); magF2=abs(F2); angF2=angle(F2); figure(4) subplot(2,2,1); plot(w/pi,magF1);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,3);plot(w/pi,angF1);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度') subplot(2,2,2); plot(w/pi,magF2);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,4);plot(w/pi,angF2);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度')
%% 折叠 F1=dtft(x,-n,w); F2=dtft(x,n,-w); magF1=abs(F1); angF1=angle(F1); magF2=abs(F2); angF2=angle(F2); figure(5) subplot(2,2,1); plot(w/pi,magF1);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,3);plot(w/pi,angF1);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度') subplot(2,2,2); plot(w/pi,magF2);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,4);plot(w/pi,angF2);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度')
%% 对称 [xf,n2]=sigfold(x,n); n3=-10:10; xe=1/2*[xf(1:10) xf(11)+x(11) x(2:11)];%奇对称分量 F1=dtft(xe,n3,w); F2=real(X); realF1=real(F1); realF2=real(F2); figure(6) subplot(2,2,1);plot(w/pi,realF1);grid xlabel('以pi为单位的频率');title('实部');ylabel('实部') subplot(2,2,2);plot(w/pi,realF2);grid xlabel('以pi为单位的频率');title('实部');ylabel('实部') xo=1/2*[-xf(1:10) -xf(11)+x(11) x(2:11)];%偶对称分量 F1=dtft(xo,n3,w); F2=j*imag(X); imagF1=imag(F1); imagF2=imag(F2); figure(7) subplot(2,2,3);plot(w/pi,imagF1);grid xlabel('以pi为单位的频率');title('虚部');ylabel('虚部') subplot(2,2,4);plot(w/pi,imagF2);grid xlabel('以pi为单位的频率');title('虚部');ylabel('虚部')
7.
w=0:pi/500:pi; H=1./(1-0.8*exp(-j*w)); magH=abs(H); angH=angle(H); figure(1) subplot(1,2,1);plot(w/pi,magH);grid xlabel('以pi为单位的频率');title('幅频特性');ylabel('幅度') subplot(1,2,2);plot(w/pi,angH);grid xlabel('以pi为单位的频率');title('相频特性');ylabel('相位')
b=1;a=[1,-0.8]; % 输入系数矩阵 n=[0:100];x=cos(0.05*pi*n); % 输入激励序列 y_zs=filter(b,a,x); % 计算响应序列 subplot(2,1,1);stem(n,x,'.'); grid % 画出激励枝干图 xlabel('n');ylabel('x(n)');title('Input sequence'); subplot(2,1,2);stem(n,y_zs,'.'); grid % 画出响应枝干图 xlabel('n');ylabel('y_zs(n)');title('Output sequence');
8.采样和重构
%% 5 a Dt=0.00005; % 模拟时间间隔：在特定精度下信号为模拟的 t=-0.005:Dt:0.005; % 模拟时刻序列 xa=exp(-1000*abs(t)); % 在特定精度下，由离散信号逼近模拟信号 fs=5000; Ts=1/fs; % 采样周期 n=-25:1:25; % 离散时间索引 x1=exp(-1000*abs(n*Ts)); % Fs=5000 样本/s:x1为采样后的离散时间序列 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的负半轴 figure(1); subplot(2,1,1);plot(t*1000,xa);hold on % 绘出xa stem(n*Ts*1000,x1,'r.'),hold off,title('时域波形') % 绘出x(jw) subplot(2,1,2);plot(w/pi,Xr);title('频域波形') % 绘出以pi归一化的数字频率对应的频域实部波形
%% 5 c fs=5000重构 Ts=0.0002; % 采样间隔 Fs=1/Ts; % 采样频率 n=-25:1:25; % 采样时间索引 nTs=n*Ts; % 序列时刻索引 Dt=0.00005; % 模拟时刻精度 t=-0.005:Dt:0.005; % 模拟时刻索引 xa= x1* ... % 内插重构 sinc(Fs*(ones(length(nTs),1)*t-nTs'*ones(1,length(t)))); figure(3); subplot(2,1,1),plot(t*1000,xa, 'k' ),hold on stem(n*Ts*1000,x1, 'r.' ),hold off ,title('重构波形(不失真)' )
9.
%% 6 (1)a Dt=0.00005; % 模拟时间间隔：在特定精度下信号为模拟的 t=-0.005:Dt:0.005; % 模拟时刻序列 xa=sin(1000*pi*t); % 在特定精度下，由离散信号逼近模拟信号 Ts=0.0001; % 采样周期 n=-25:1:25; % 离散时间索引 x1=sin(1000*pi*abs(n*Ts)); % Fs=5000 样本/s:x1为采样后的离散时间序列 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的负半轴 figure(); subplot(2,2,1);plot(t*1000,xa);hold on % 绘出xa stem(n*Ts*1000,x1,'r.'),hold off,title('时域波形') % 绘出x(jw) subplot(2,2,2);plot(w/pi,Xr);title('频域波形') % 绘出以pi归一化的数字频率对应的频域实部波形 magXr=abs(Xr); angXr=angle(Xr); subplot(2,2,3); plot(w/pi,magXr);grid xlabel('以pi为单位的频率');title('幅度部分');ylabel('幅度') subplot(2,2,4);plot(w/pi,angXr);grid xlabel('以pi为单位的频率');title('相角部分');ylabel('弧度')
%% 6 (2)a Dt=0.00005; % 模拟时间间隔：在特定精度下信号为模拟的 t=0:Dt:1; % 模拟时刻序列 xa=sin(20*pi*t+pi/4); % 在特定精度下，由离散信号逼近模拟信号 Ts=0.05; % 采样周期 n=0:1:20; % 离散时间索引 x1=sin(20*pi*abs(n*Ts)+pi/4); % Fs=5000 样本/s:x1为采样后的离散时间序列 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的负半轴 figure(); plot(t,xa);hold on % 绘出xa stem(n*Ts,x1,'o'),hold off,title('时域波形') % 绘出x(jw)
%% 6 (2)b Ts=0.05; % 采样间隔 Fs=1/Ts; % 采样频率 n=0:1:20; % 采样时间索引 nTs=n*Ts; % 序列时刻索引 Dt=0.001; % 模拟时刻精度 t=0:Dt:1; % 模拟时刻索引 xa= x1* ... % 内插重构 sinc(Fs*(ones(length(nTs),1)*t-nTs'*ones(1,length(t)))); figure(); plot(t,xa, 'k' ),hold on stem(n*Ts,x1, 'r.' ),hold off ,title('重构波形(不失真)' )
