%經典譜估計直接法實現FSK信號的解調
%中頻fc為15MHz，采樣率fs為50MHz
%一個脈沖內5bit信息
%fs=n/tb，tb為碼元寬度,n為一個碼元內的采樣點數
clear all;

fs=5.0e+007;
ts=1/fs;
tb=5e-006;%碼元寬度5us
sample=tb/ts;%一個碼元內采樣250個點
fc=1.5e+007;
f1=fc+8/(tb);%1.66e+007,對應1
f2=fc-8/(tb);%1.34e+007，對應0
B=5.0e+006;


%對該輸入信號FSK調制
Xin=[1,1,1,1,1,1,0,1,1,1];%randn(5,1)>0;
tt=(0:tb/sample:tb-(tb/sample));
t=[tt;tt+tb;tt+2*tb;tt+3*tb;tt+4*tb;tt+5*tb;tt+6*tb;tt+7*tb;tt+8*tb;tt+9*tb];
y=zeros(10,length(tt));
for i=1:10
y(i,:)=Xin(i)*cos(2*pi*f1*t(i,:))+~Xin(i)*cos(2*pi*f2*t(i,:));
end
t=reshape(t',length(tt)*10,1);
y=reshape(y',length(tt)*10,1);


%變換到基帶
p1=cos(2*pi*fc*t-2*B*pi*t/2);
z=y.*p1;


%窗設計低通FIR濾波器
%假設信號帶寬為3.5MHz
fs=3.22e+007;
fn=1.0e+007;%截止頻率
Wn=fn/fs;
N=127;%低通FIR濾波器的階數
K=FIR1(N,Wn);
H=zeros(1,2500);
Y=conv(z,K);
n=250;
%截取后320個點
for j=1:2500
    H(j)=Y(j+63);
end
%長度為1250

%解調FSK信號
%譜估計：直接法
k1=floor(250*(f1-fc+B/2)/fs)
k2=floor(250*(f2-fc+B/2)/fs)
for i=1:10
    for j=1:n
        C(i,j)=H(j+(i-1)*n);
    end
end

OUT=zeros(1,10);
for i=1:10
    A=fft(C(i,:));
    for j=1:n
      Pxx(j)=(abs(A(j))^2)/n;
    end
    if(Pxx(k1)>Pxx(k2))
               OUT(i)=1;
    else
               OUT(i)=0;
    end
end
OUT