ddc_simu.m

用MATLAB实现软件无线电中DDC的仿真

%仿真时钟抖动对下变频的影响  用76.8采96中频信号
%几个图分别是
%s1时域采样，下变频，采样时间t=n*t0加入了时间抖动，抖动用rand函数产生
%s2下变频后，下变频，
%s3 60倍抽取，
%s4 FFT，
%s5 FFT移到中间，
%s9 76.8m时钟加入抖动，不过这个加抖动的方法不是很有道理，s9的产生可再考虑
%满足N*fs/(fre-96)=num,N为整数时,不用加窗
%满足N为整数,且num与N互质,这种情况下SNR最大(比其他情况大1dB左右)
%加权对SNR的损失最大4-5dB
%未考虑采样位数的影响
%
clear all;
clc;
close all;
format long;
fre=96.02;
fre2=76.8;
nco=19.2;
t0=1/76.8;
fhs=76.8;
fs=1.28;
rati=5000; %determine clock jitter 
num=128*60;
  jitter=(rand(1,num)-0.5)/rati;
  for i=1:num
  %    jitter(i)=0;
  end
for n=1:num
  s1(n)=200*exp(j*2*pi*fre*(n+jitter(n)-1)*t0);
  s2(n)=s1(n)*exp(-j*2*pi*nco*(n-1)*t0);
 % s9(n)=200*exp(j*2*pi*fre2*(n+jitter(n)-1)*t0);
end
for n=1:num/60
 s3(n)=s2(60*n);
end
jitter_rms=0;
for n=1:num
jitter_rms=jitter(n)*jitter(n)+jitter_rms;
end
jitter_rms=sqrt(jitter_rms/num)*1000*1000*t0

figure(1);
plot(1:num,real((s1(1:num))));
title('中频采样信号');
figure(2);

plot(1:num,real((s2(1:num))));
title('下变频后信号');

figure(3);

plot(1:num/60,real((s3(1:num/60))));

title('60倍抽取后信号');

fftp=num/60;

w = hamming(num/60);                                % hamming weight
w=chebwin(num/60,80);
%s3 = s3.*w'; 			 % comment this instrument if you don't want hamming weight

s4=fft(s3,fftp);

[pp,ii]=max(abs(s4));
pp=0;
for i=ii-0:ii+3
    pp=abs(s4(i))*abs(s4(i))+pp;
end

qq=0;
for i=1:fftp
    qq=abs(s4(i))*abs(s4(i))+qq;
end
pp
qq
snr = 10*log10(pp/(qq-pp));
for n=(fftp/2+1):fftp
 s5(n-fftp/2)=s4(n);
end
for n=1:fftp/2
 s5(fftp/2+n)=s4(n);
end
figure(4);

stem(1:fftp,(abs(s5)));
title('采样信号的FFT');

figure(5);
f = fs*(0:fftp/2-1)/fftp;
ffff=t0/rati*1000*1000;
snr2=-20*log10(2*pi*fre*jitter_rms*1e-6);       %using adi formula
xxxxxx=sprintf('加入jitter后的采样信号的频谱 jitter RMS is %fps SNR is %fdB ADI result is %fdB',jitter_rms,snr,snr2);
stem(f,20*log10(abs(s4(1:fftp/2))),'.');
xlabel('frequency (MHz)')

title(xxxxxx);


fre2=76.8;
figure(6);
t0=0.02/fre2

num=500;
   rati=100;%jitter 的单位  =1/rati  ns
  jitter=(rand(1,num)-0.5)/rati;
  for i=1:num
    % jitter(i)=0;
  end
for n=1:num
  s9(n)=2000*exp(j*2*pi*fre2*(n+jitter(n)-1)*t0);
end

w=chebwin(num,100);
w=hamming(num);

%s9 = s9.*w'; 									      % comment this instrument if you don't want hamming weight

jitter_rms=0;
for n=1:num
jitter_rms=jitter(n)*jitter(n)+jitter_rms;
end
jitter_rms=sqrt(jitter_rms/num)*1000*1000*t0;         %jitter_rms受采样间隔影响,所以即使rati相同,jitter_rms也可能不同
xxxxxx=sprintf('加入jitter后的时钟信号的频谱 jitter RMS is %fps',jitter_rms);

plot(1:num,real(s9(1:num)));
title('加入jitter后的时钟信号');

fs=1/t0;
fftp=num;
figure(7);
f = fs*(0:fftp/2-1)/fftp;
ffts9=fft(s9,fftp);
stem(f,20*log10(abs(ffts9(1:fftp/2))),'.');
title(xxxxxx);