D=8;fs=800*10^3;ts=1/fs;
t=0:ts:0.001;
s=(2.*cos(2*pi*15*10^3*t)+2.*cos(2*pi*65*10^3*t)+3.*cos(2*pi*115*10^3*t)+4.*cos(2*pi*165*10^3*t)+3.*cos(2*pi*215*10^3*t)+2.*cos(2*pi*265*10^3*t)+2.*cos(2*pi*315*10^3*t)+3.*cos(2*pi*365*10^3*t)).*cos(2*pi*4*10^6*t);%經過采樣的射頻模擬信號

plot(t,s);%時域圖象
pause

%功率譜
fs0=10*10^6;ts0=1/fs0;
t0=0:ts0:0.001;
s0=(2.*cos(2*pi*15*10^3*t0)+2.*cos(2*pi*65*10^3*t0)+3.*cos(2*pi*115*10^3*t0)+4.*cos(2*pi*165*10^3*t0)+3.*cos(2*pi*215*10^3*t0)+2.*cos(2*pi*265*10^3*t0)+2.*cos(2*pi*315*10^3*t0)+3.*cos(2*pi*365*10^3*t0)).*cos(2*pi*4*10^6*t0);%經過采樣的射頻模擬信號

signal=[t0' s0'];
a=[s0'];
y=fft(s0,512);
p=y.*conj(y)/512;
f=fs0*(0:256)/512;
plot(f,p(1:257));
pause

%低通濾波器設計
rp = 3;         
rs = 40;         
fs = 800000;      
f = [370000 400000];   
a = [1 0];      
dev = [(10^(rp/20)-1)/(10^(rp/20)+1)  10^(-rs/20)]; 
[n,fo,ao,w] = remezord(f,a,dev,fs);
b = remez(n,fo,ao,w);
freqz(b,1,1024,fs);
title('Lowpass Filter Designed to Specifications');
pause

%信號通過低通濾波器
sn=conv(s,b);

%信號通過低通濾波器后的功率譜
fs00=800*10^3;
y0=fft(sn,512);
p00=y0.*conj(y0)/512;
f00=fs00*(0:256)/512;
plot(f00,p00(1:257));
pause


%原型低通濾波器設計
rp = 2;         
rs = 62;         
fs = 800000;       
f = [20000 25000];   
a = [1 0];       

dev = [(10^(rp/20)-1)/(10^(rp/20)+1)  10^(-rs/20)]; 
[n,fo,ao,w] = remezord(f,a,dev,fs);
b = remez(n,fo,ao,w);
freqz(b,1);
title('Lowpass Filter Designed to Specifications');
pause



%信號的信道化設計
for k=1:D
    for n=8:830
        szn(k,n-7)=sn(n-(k-1));
    end
end     


for k=1:D
   for m=1:103
    szm(k,m)=szn(k,m*8-7);
    szxm(k,m)=szm(k,m)*((-1)^m)*exp(j*(pi/2)*m);
    
    end    
end 


%每一支路的濾波器
for r=1:36
    for k=1:D
        h(k,r)=b((r-1)*D+k);
    end
end    

%濾波
for k=1:D
    for r=1:36
        hk(r)=h(k,r);
    end
    
    for m=1:103
        szxmk(m)=szxm(k,m);
    end
    
     szxhmk=conv(szxmk,hk);
    for m=1:138
        szxhm(k,m)=szxhmk(m);    
    end
end    

for k=1:D
    for m=1:69
        szxh2m(k,m)=szxhm(k,2*m-1);
        szxh2xm(k,m)=szxh2m(k,m)*(-1)^(k-1)*exp(-j*(pi/16)*(k-1));
    end
end

%DFT
for m=1:69
    for k=1:D
        xk(k)=szxh2xm(k,m);
    end
        xfft=fft(xk);
   
    for k=1:D;
    y1111(k,m)=xfft(k);
    end
end

%信道1內的功率譜
fs11=100*10^3;
for m=1:69
    y1(m)=y1111(1,m);
end
y11=fft(y1,512);
p11=y11.*conj(y11)/512;
f11=fs11*(0:256)/512;
plot(f11,p11(1:257));
pause

%信道2內的功率譜
for m=1:69
    y1(m)=y1111(2,m);
end
y11=fft(y1,512);
p11=y11.*conj(y11)/512;
f11=fs11*(0:256)/512;
plot(f11,p11(1:257));
pause

%信道3內的功率譜
for m=1:69
    y1(m)=y1111(3,m);
end
y11=fft(y1,512);
p11=y11.*conj(y11)/512;
f11=fs11*(0:256)/512;
plot(f11,p11(1:257));
pause

%信道4內的功率譜
for m=1:69
    y1(m)=y1111(4,m);
end
y11=fft(y1,512);
p11=y11.*conj(y11)/512;
f11=fs11*(0:256)/512;
plot(f11,p11(1:257));
pause

%信道5內的功率譜
for m=1:69
    y1(m)=y1111(5,m);
end
y11=fft(y1,512);
p11=y11.*conj(y11)/512;
f11=fs11*(0:256)/512;
plot(f11,p11(1:257));
pause

%信道6內的功率譜
for m=1:69
    y1(m)=y1111(6,m);
end
y11=fft(y1,512);
p11=y11.*conj(y11)/512;
f11=fs11*(0:256)/512;
plot(f11,p11(1:257));
pause

%信道7內的功率譜
for m=1:69
    y1(m)=y1111(7,m);
end
y11=fft(y1,512);
p11=y11.*conj(y11)/512;
f11=fs11*(0:256)/512;
plot(f11,p11(1:257));
pause

%信道8內的功率譜
for m=1:69
    y1(m)=y1111(8,m);
end
y11=fft(y1,512);
p11=y11.*conj(y11)/512;
f11=fs11*(0:256)/512;
plot(f11,p11(1:257));
pause