fm=1;
fc=10;
t=0:0.01:5;
mt=sin(2*pi*fm*t);
ct=cos(2*pi*fc*t);
am=(1+0.8*mt).*ct;
pm=cos(2*pi*fm*t);
fm=cos(2*pi*fc*t+2*pi*0.8*pm);
plot(t,ct+9,t,mt+6,t,am+3,t,fm);
axis([0 4 -2 11])

