clear all;
close all;
kf=5;
fc=2.5;
T=5;dt=0.001;
t=0:dt:T;
%信源
fm=1;fs=32;fsc=40;
%mt=cos(2*pi*fm*t)+1.5*sin(2*pi*0.3*fm*t);%信源信號
mt=cos(2*pi*fm*t);%信源信號
%fm調制
A=sqrt(2);
%mti=1/2/pi/fm*sin(2*pi*fm*t)-3/4/pi/0.3/fm*cos(2*pi*0.3*fm*t);
mti=1/2/pi/fm*sin(2*pi*fm*t);%mt的積分函數
st=A*cos(2*pi*fc*t+2*pi*kf*mti);
figure(2)
subplot(311)
plot(t,st);hold on;
plot(t,mt,'r--');
xlabel('t');ylabel('調頻信號');
subplot(312)
[f,sf]=T2F(t,st);
plot(f,abs(sf));
axis([-25 25 0 3])
xlabel('f');ylabel('調頻信號幅度譜')
%fm解調
for k=1:length(st)-1
    rt(k)=(st(k+1)-st(k))/dt;
end
rt(length(st))=0;
subplot(313)
plot(t,rt);hold on;
plot(t,A*2*pi*kf*mt+A*2*pi*fc,'r--');
xlabel('t');ylabel('調頻信號微分后包絡')

N=3072;    %數據長度
M=64;        %頻率平滑長度
figure(3);
K=round(2.5*N*fc/fsc);          %頻率軸
alfa=round(2.5*N*fc/fsc);
G1=pxg(st,N,M,K,alfa,fs);
p=(-K:32:K)*fsc/N;
q=(-alfa:32:alfa)*2*fsc/N;
mesh(p,q,abs(G1));