离散系统的变换域分析【数字信号处理四】

一、求系统H(z)的零、极点、幅频响应和相位响应

num=[0.0528,0.0797,0.1295,0.1295,0.797,0.0528];
den=[1,-1.8107,2.4947,-1.8801,0.9537,-0.2336];
k=2048;
w=0:pi/k:pi;
h=freqz(num,den,w);
figure(1);
zplane(num,den); 
title("零极点图")
figure(2);
subplot(2,1,1);
plot(w/pi,abs(h));
title('幅度谱'); 
xlabel('\omega/\pi');ylabel('幅值');
subplot(2,1,2);
plot(w/pi,angle(h));
title('相位谱'); 
xlabel('\omega/\pi');ylabel('弧度');

二、用Matlab验证DFT运算的对称性质

N=32;
n=0:N-1;
nn=mod(-n,N)+1;
xn=cos(2*pi*6*n/N)+1i*sin(2*pi*6.8*n/N);
% x[n]  -->  X[k]
Xk=fft(xn);
subplot(2,7,1)
stem(n,abs(Xk));
title("x[n]  -->  X[k]")

% x*[n]  -->  X*[<-k>N]
xn1=conj(xn);
Xk1=fft(xn1);
subplot(2,7,2)
stem(n,abs(Xk1));
title("x*[n]  -->  X*[<-k>N]");
subplot(2,7,2+7);
stem(n,abs(conj(Xk(nn))));
title("直接求X*[<-k>N]");

% x*[<-n>N]  -->  X*[k]
xn2=conj(xn(nn));
Xk2=fft(xn2);
subplot(2,7,3);
stem(n,Xk2);
title("x*[<-n>N]  -->  X*[k]");
subplot(2,7,3+7);
stem(n,conj(Xk));
title("直接求X*[k]");

% xre[n]  -->  Xcs[k]=1/2{X[k]+X*[<-k>N]}
xn3=real(xn);
Xk3=fft(xn3);
subplot(2,7,4)
stem(n,abs(Xk3));
title("xre[n]  -->  Xcs[k]=1/2(X[k]+X*[<-k>N])");
subplot(2,7,4+7);
stem(n,1/2*(Xk+conj(Xk(nn))));
title("直接求Xcs[k]");

% xim[n]  -->  Xca[k]=1/2{X[k]-X*[<-k>N]}
xn4=1i*imag(xn);
Xk4=fft(xn4);
subplot(2,7,5)
stem(n,Xk4);
title("xim[n]  -->  Xca[k]=1/2(X[k]-X*[<-k>N])");
subplot(2,7,5+7);
stem(n,1/2*(Xk-conj(Xk(nn))));
title("直接求Xca[k]");

% xcs[n]  -->  Xre[k]
xncs=1/2*(xn+conj(xn(nn)));
Xk5=fft(xncs);
subplot(2,7,6)
stem(n,Xk5);
title("xim[n]  -->  Xre[k]");
subplot(2,7,6+7)
stem(n,real(Xk));
title("直接求Xre[k]");

% xca[n]  -->  jXim[k]
xnca=1/2*(xn-conj(xn(nn)));
Xk6=fft(xnca);
subplot(2,7,7)
stem(n,Xk6);
title("xca[n]  -->  jXim[k]");
subplot(2,7,7+7)
stem(n,1i*imag(Xk));
title("直接求jXim[k]");

第一列为求x(n)变化后对应的fft，第二列为直接求fft。
结果一致，证明了DFT运算的对称性质。

三、产生数字线性调频信号，分析时域波形和频谱特性

fre=5*10.^7;
T=1/fre;
t=0.1*10.^-3;
fre_start=9.5*10.^5;
fre_end=1.05*10.^6;
k=2*(fre_end-fre_start)/t;

n=0:T:t;
len=length(n);
x=cos(2*pi*(fre_start+1/2*k*n).*n);
subplot(2,1,1);
plot(n,x);
xk=abs(fftshift(fft(x)));
subplot(2,1,2);
f=(-len/2:len/2-1)*fre/len;
stem(f,xk);
axis([0,2*10.^6,0,800]);
title('频谱图');xlabel('f/Hz');

时域上基本看不出调频后有什么变化。
但观察频谱可以发现，通过调频，频谱上有了对应的分量，范围在1Ghz左右，符合理论计算结果。

四、设计简单的OFDM系统，并验证循环前缀的作用

M=8;N=3;
xn=round(rand(M,1));
hn=[1,0.1,0.01];

xnk=ifft(xn);
xncp=[xnk(end-(N-2):end,1);xnk];

yncp=conv(xncp,hn);
Yk=yncp(N:N+M-1).';

hncp=[hn,zeros(1,M-N)];
Hk=fft(hncp);

for i=1:M
    Xk(i)=Yk(i)/Hk(i);
end
xnifft=abs(fft(Xk));
err=0.5;
for i=1:M
    if xnifft(i)<err
        xn_output(i)=0;
    else
        xn_output(i)=1;
    end
end
k=0:M-1;
subplot(3,1,1);
stem(k,xn);
title('随机生成数列xn');
subplot(3,1,2);
stem(k,xnifft);
title('接收后的yn序列');
subplot(3,1,3);
stem(k,xn_output);
title('设置阈值处理误差后恢复的yn序列');

随机序列xn长度为M，信道冲激响应hn长度为N。
首先我们对随机序列xn进行IDFT变换，然后通过N-1个样本添加循环前缀将xn延长为M+N-1的序列xncp。
将xncp与hn卷积得到yncp，取中间M个样本得到yk。
将hn补零至M点求FFT得到hk
Xk等于yk / hk
对Xk求FFT即得到yn
此时的yn有误差，通过设置阈值，判决处理后恢复了xn。