%數字下變頻全過程仿真程序 寬帶

%fir1
%信號帶寬為30MHz
%輸入采樣率為96MHz
%輸出采樣率為32MHz，抽取3

clear all;
close all;
fsamp=102.4e6;
Ts=1/fsamp;
fc=32.4e6;
band=30e6;
Tp=40e-6;
N=Tp*fsamp;
n=0:N-1;
u=band/Tp; %u為帶寬除以時寬
t=-Tp/2:Tp/N:Tp/2-Tp/N;
f0=70e6;


xs=cos(2*pi*(f0*t+0.5*u*t.^2));         %信號表達式
S0=fft(xs,N);
S1=abs(S0);
S2=(S1);         %信號頻譜
ll=length(S2);
S2=awgn(S2,10);
f=0:fsamp/N:fsamp-fsamp/N;
figure(1);
plot(f,20*log10(S2/max(S2)));title('輸入LPM（30M 40us）信號頻譜');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');
grid on

%----------generate NCO---------%
for t=1:N
    t1=(t-1)*Ts;
    ncoi_c(t)=cos(2*pi*fc*t1);
end

for t=1:N
    t1=(t-1)*Ts;
    ncoq_c(t)=sin(2*pi*fc*t1);
end
ncoi=awgn(ncoi_c,80);
ncoq=awgn(ncoq_c,80);

f=0:fsamp/N:fsamp-fsamp/N;
u1=abs(fft(ncoi));
u2=abs(fft(ncoq));
figure(2);plot(f,20*log10(u1/max(u1)));title('NCO（I 路)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');
grid on
figure(3);plot(f,20*log10(u2/max(u2)));title('NCO（Q 路)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');
grid on


for n=1:1:N
    ysi(n)=xs(n)*ncoi(n);
    ysq(n)=xs(n)*ncoq(n);
end


u1=abs(fft(ysi));
u2=abs(fft(ysq));
f=0:fsamp/N:fsamp-fsamp/N;
figure(4);plot(f,20*log10(u1/max(u1)));title('混頻輸出（I 路)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');
grid on
figure(5);plot(f,20*log10(u2/max(u2)));title('混頻輸出（Q 路)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');
grid on
% u=u1+u2;
u=abs(fft(ysi-j*ysq));
figure(6);plot(f,20*log10(u/max(u)));title('混頻輸出（I+JQ 路)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');
grid on

%設計fir1濾波器
fs=fsamp;
band=30e6;      %信號帶寬
gdd=7.5e6;    %過渡帶寬 
wp=(band/2)*2*pi/fs;ws=(band/2+gdd)*2*pi/fs;
tr_width=gdd*2*pi/fs;   %過渡帶寬
rp=0.5;  %通帶波紋小于0.01dB
rs=90;  %最小阻帶衰減80dB
f=[band/2 band/2+gdd]; %理想低通濾波器的截止頻率
a=[1 0];  %Desired amplitudes
dev=[(10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)];
[NN,fo,ao,w]=remezord(f,a,dev,fs);    %NN為濾波器長度
NN=41;
fira=remez(NN,fo,ao,w);
figure(7);
freqz(fira,1,1024,fs);
title('抗混疊FIR濾波器');%[hf,f1]=freqz(fira,1,NN,fs);

% figure(7);
% f=f1;
% adg=1*20*log10(abs(hf));
% plot(f,adg);grid;title('FIR1 Magnitude Response(dB)');
% xlabel('frequency (Hz)');
% ylabel('Magnitude (dB)');
N_FIR1=NN

ysi1=conv(fira,ysi);
ysq1=conv(fira,ysq);

Nu=length(ysi1);
f=0:fs/Nu:fs-fs/Nu;
figure(8);
adg=abs(fft(ysi1));
plot(f,20*log10(adg/max(adg)));grid;title('FIR1  output  Magnitude Response(dB) (I)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');
figure(9);
ysi_IQ1=ysi1-j*ysq1;
adg=abs(fft(ysi_IQ1));
plot(f,20*log10(adg/max(adg)));grid;title('FIR輸出信號頻譜(I+JQ)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');

ysi11(1:ceil(Nu/2))=ysi1(1:2:Nu);
ysq11(1:ceil(Nu/2))=ysq1(1:2:Nu);

figure(10);
Nu=length(ysi11);
adg=abs(fft(ysi11));
f=0:fs/(2*Nu):fs/2-fs/(2*Nu);
plot(f,20*log10(adg/max(adg)));grid;title('dec3 Magnitude Response(dB) (I)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');

ysout=ysi11-j*ysq11;

figure(12);
Nu=length(ysout);
adg=abs(fft(ysout));
f=0:fs/(2*Nu):fs/2-fs/(2*Nu);
plot(f,20*log10(adg/max(adg)));grid;title('DDC輸出信號頻譜(dB)(I+JQ)');
xlabel('frequency (Hz)');
ylabel('Magnitude (dB)');