数字信号处理第三次实验代码及波形图

数字信号处理第三次实验
一、实验目的：
1、掌握离散时间系统的DFT的MATLAB实现；
2、熟悉DTFT和DFT之间的关系。

3、了解信号不同变形的DFT与原信号DFT之间的关系
二、实验内容：
1．开发一个时域圆周移位的matlab函数cirshift.m并测试：设x(n)=[9 8 7 6 5 4 3],求x((n-4))8R8(n)以及x((n+5))8R8(n)。

function y=cirshift(x,m,N)
%长度为N的x序列（时域）作m点圆周位移
%-------------------------------------
%[y]=cirshift(x,m,N)
%y=包含圆周位移的输出序列
%x=长度<=N的输入序列
%m=移位样点数
%N=圆周缓冲器长度
%方法：y(n)=x((n-m)mod N)
%check for length of x
if length(x)>N
error('N必须>=x的长度')
end
x=[x zeros(1,N-length(x))];%将x补零到长度为N
n=[0:1:N-1];
n=mod(n-m,N);
y=x(n+1);
%主函数.m
clc;clear all;
x=[9 8 7 6 5 4 3];
y1=cirshift(x,4,8);
y2=cirshift(x,-5,8);
2. 开发一个圆周翻转的matlab函数cirflip.m并测试：设x(n)=[9 8 7 6 5 4 3 2 1],求x(n)的16点圆周翻转x((-n))16R16(n)。

function y=cirfilp(x,N)
nx=[0:1:N-1];
y=x(mod(-nx,N)+1);
%主函数.m
clc;clear all;
x=[9 8 7 6 5 4 3 2 1];
N=16;
x=[x zeros(1,N-length(x))];
nx=0:N-1;
y=cirfilp(x,16);
subplot(121),stem([0:N-1],x);title('原序列');
xlabel('n');ylabel('x(n)');grid;
subplot(122),stem([0:N-1],y);title('圆周翻褶序列');
xlabel('n');ylabel('x((n))16r16(n)');grid;
3. 开发一个分解圆周奇偶对称序列的matlab函数cirevod.m并测试：设x(n)=[9 8 7 6 5 4 3 2 1],将x(n)分解圆周奇偶对称序列
x ep(n) 、x op(n)。

function [xep,xop]=circevod(x)
%将实序列分解为圆周偶和圆周奇两部分
%------------------------------------
%[xep,xop]=circevod(x)
if any(imag(x)~=0)
error('x非实序列')
end
N=length(x);
n=0:N-1;
xep=0.5*(x+x(mod(-n,N)+1));
xop=0.5*(x-x(mod(-n,N)+1));
%主函数.m
x=[9 8 7 6 5 4 3 2 1];
[xep,xop]=circevod(x);
4. 开发一个实现DFT矩阵的函数wN=dftmtx(N)并测试N=4，8，16的DFT矩阵，在此基础上寻求IDFT矩阵wNI的编程方法并测试N=4，8，16的IDFT矩阵
%这是matlab自带的函数命令行输入type+dftmtx可查看源代码。

function D = dftmtx(n)
n = signal.internal.sigcasttofloat(n,'double','dftmtx','N',...
'allownumeric');
D = fft(eye(n));
命令行输入：
wN=dftmtx(N);
wNI=conj(dftmtx(N))/N
DFT
N=4
N=8
N=16
IDFT：N=4
N=8
N=16
5.自编工具函数XK=dft(xn,N)实现DFT以及IDFT，并测试x(n)=[9 8 7 6 5 4 3 2 1]的DFT。

function[Xk]=dft(xn,N)
n=[0:1:N-1];
k=n;
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(nk);
Xk=xn*WNnk;
%主函数
xn=[9 8 7 6 5 4 3 2 1]
y=dft(x,9);
function[Xk]=idft(xn,N)
n=[0:1:N-1];
k=n;
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(-nk);
Xk=(xn*WNnk)/N;
6. 开发一个圆周卷积计算的matlab函数y=circonvt(x1,x2,N),并测试x1(n)=[1 2 3 4 9], x2=[7 6 5 4 3 2 1],求8点圆周卷积。

function y=circonvt(x1,x2,N)
if length(x1)>N
error('N必须>=x1的长度')
end
if length(x2)>N
error('N必须>=x2的长度')
end
x1=[x1 zeros(1,N-length(x1))]; x2=[x2 zeros(1,N-length(x2))]; m=[0:1:N-1]; x2=x2(mod(-m,N)+1); H=zeros(N,N);
for n=1:1:N
H(n,:)=cirshift(x2,n-1,N); end
y=x1*H';
%test.6
x1=[1 2 3 4 9];
x2=[7 6 5 4 3 2 1];
y=circonvt(x1,x2,8);
7.教材P212 习题3.36.
余弦频率是正弦频率的两倍，故正弦周期是预先周期的两倍。

T=0.5s 抽样间隔T=0.01s，故N=50；
clear all;clc;
N=50;
n=[0:1:N-1];
k=n;
xn=2*sin(0.04*pi*n)+5*cos(0.08*pi*n);
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(nk);
Xk=xn*WNnk;
magX=abs(Xk);
angX=angle(Xk);
subplot(121);stem(k,magX,'.');grid;
xlabel('k');title('幅度谱曲线');ylabel('X(k)');
subplot(122);stem(k,angX,'.');grid;
xlabel('k');title('相位谱曲线');ylabel('弧度');
因为抽样频率不变，要增大N才能减少谱泄露，故N=99；clear all;clc;
N=99;
n=[0:1:N-1];
k=n;
xn=2*sin(0.04*pi*n)+5*cos(0.08*pi*n);
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(nk);
Xk=xn*WNnk;
magX=abs(Xk);
angX=angle(Xk);
subplot(121);stem(k,magX,'.');grid;
xlabel('k');title('幅度谱曲线');ylabel('X(k)'); subplot(122);stem(k,angX,'.');grid;
xlabel('k');title('相位谱曲线');ylabel('弧度');。