%數字下變頻仿真程序
%數字下變頻的主要作用有兩點：

%一、對A/D采樣后的高頻/中頻信號序列進行頻譜搬移（通過與數控振蕩器產生的數字本振信號序列進行相乘下變頻到基頻）。
%二、對基頻上高采樣率的信號序列進行抽取，多采樣率變換，降低數字信號序列密度。
??????????????????????????
%%%%%實際的數字下變頻在對高頻/中頻信號序列進行A/D采樣之前為了防止發生頻率混疊，要進行預濾波處理。
%%%%在對基頻上高采樣率的信號序列進行抽取之前，要通過CIC濾波，HB濾波，FIR濾波，以防止抽取時發生頻譜混疊。
?????????????????????????
????????????????????????
clear?all;
close?all;
fsamp=96e6;????????%fsamp=96MHz?輸入采樣頻率為96MHz
Ts=1/fsamp;????????%Ts為fsamp的倒數?即輸入采樣間隔Ts
band=30e6;?????????%預設的采樣帶寬為30MHz
Tp=60e-6;??????????%預設的采樣時間周期Tp為60us
N=Tp*fsamp;????????%N為輸入采樣頻率與采樣時間周期之積。表示在采樣時間周期Tp內，以fsamp的采樣率采樣可以得到的采樣點數???N?=?5760
u=band/Tp;?????????%u為帶寬除以時寬。表示在單位時間間隔內的頻帶寬度。也即這30M的帶寬分布在Tp=60us的時間周期上，單位時間的頻帶寬度
t=-Tp/2:Tp/N:Tp/2-Tp/N;???????%t取點從-Tp/2開始以Tp/N為步進值增加到Tp/2-Tp/N。

f0=70e6;???????????%輸入的已調頻信號載波頻率為70MHz
xs=cos（2*pi*（f0*t+0.5*u*t.^2））;????????%輸入的已調頻信號經fsamp=96MHz帶通采樣后的輸出。相當于A/D轉換后的數字信號序列
S0=fft（xsN）;??????%S0是對A/D轉換后的數字信號序列進行N點fft的結果；
S1=abs（S0）;????????%S1是對s0求模的結果；
S2=（S1）;
S2=awgn（S210）;????%在S2的頻譜中加入10db的高斯白噪聲
f=0:fsamp/N:fsamp-fsamp/N;???????%f的取點由0開始以fsamp/N為步進值直到fsamp-fsamp/N結束
figure（1）;?????????%畫出載頻為70MHz有用信號帶寬為30MHz的帶通信號序列頻譜
plot（f/1e620*log10（S2/max（S2）））;???????????%橫軸以MHz為單位，縱軸以dB形式，其中S2/max（S2）表示輸出該帶通信號序列相對幅度大小，對它取對數后的結果就是dB的形式了。
title（‘載頻為70MHz有用信號帶寬為30MHz的帶通信號序列頻譜‘）;
xlabel（‘frequency（MHz）‘）;???????????????????%從頻譜圖中可以看出信號序列頻譜有（11MHz，41MHz）和（55MHz，85MHz）兩部分。
ylabel（‘Magnitude（dB）‘）;????????????????????%這是因為經fsamp=94MHz帶通采樣時，在fsamp/2=48MHz處發生了頻譜折疊，原來信號序列頻譜（55M，85M）折疊到（11MHz，41MHz）了。
grid?on?????????????????????????????????????%這兩部分頻譜形狀一致，沒有發生頻譜混疊。這里的帶通采樣速率fsamp=96MHz是通過計算得出來的。
????????????????????????????????????????????%具體計算式如下：fsamp>=4f0/（2n+1）?且fsamp>=2B。這里f0指信號序列中心頻率，B指信號序列帶寬?n取正整數。
?
%NCO數控振蕩器模塊%
for?t=1:N
????t1=（t-1）*Ts;
????ncoi_c（t）=cos（2*pi*f0*t1）;?????????????%產生頻率為f0的cos數控本振（I路），這里產生數控本振的時間間隔與A/D采樣間隔相同，便于序列后續相乘。
end
?
for?t=1:N
????t1=（t-1）*Ts;
????ncoq_c（t）=sin（2*pi*f0*t1）;????????????%產生頻率為f0的sin數控本振（Q路），這里產生數控本振的時間間隔與A/D采樣間隔相同，便于序列后續相乘。
end
ncoi=awgn（ncoi_c80）;?????????%對產生的I路本振序列加入80dB的高斯白噪聲
ncoq=awgn（ncoq_c80）;?????????%對產生的Q路本振序列加入80dB的高斯白噪聲
?
f=0:fsamp/N:fsamp-fsamp/N;????%f取點從0開始以fsamp/N為步進值直到fsamp-fsamp/N結束
u1=abs（fft（ncoi））;????????????%對加入高斯白噪聲的I路本振信號序列進行FFT后取模
u2=abs（fft（ncoq））;????????????%對加入高斯白噪聲的Q路本振信號序列進行FFT后取模
figure（2）;????????????????????%畫出加入高斯白噪聲后的數控振蕩器I路信號頻譜
plot（f/1e620*log10（u1/max（u1）））;%橫軸以MHz為單位，縱軸是dB形式
title（‘加入高斯白噪聲的數控振蕩器I路信號頻譜‘）;
xlabel（‘frequency（MHz）‘）;
ylabel（‘Magnitude（dB）‘）;
grid?on
figure（3）;?????????????????????%畫出加入高斯白噪聲后的數控振蕩器Q路信號頻譜
plot（f/1e620*log10（u2/max（u2）））;%橫軸以MHz為單位，縱軸是dB形式
title（‘加入高斯白噪聲的數控振蕩器Q路信號頻譜‘）;
xlabel（‘frequency（MHz）‘）;
ylabel（‘Magnitude（dB）‘）;
grid?on
?
for?n=1:1:N
????ysi（n）=xs（n）*ncoi（n）;%A/D帶通采樣后信號序列與數字本振I路信號序列混頻相乘（下變頻）
????ysq（n）=xs（n）*ncoq（n）;%A/D帶通采樣后信號序列與數字本振Q路信號序列混頻相乘（下變頻）
end
?
u1=abs（fft（ysi））;%對I路下變頻后序列進行FFT，取模后結果送到u1
u2=abs（fft（ysq））;%對Q路下變頻后序列進行FFT，取模后結果送到u2
f=0:fsamp/N:fsamp-fsamp/N;%f取點從0開始，以fsamp/N為步進值直到fsamp-fsamp