Basically, your code looks fine - the problem was that the FFT resolution was
not sufficient to resolve the demodulated signal with higher sample rates (which
you have set at 10 times the carrier). See my modified code below. I set L =
100*fc. I also added a windowing function prior to the FFT to reduce edge
effects. Also, make sure you restrict your parameters so as to avoid violating
Nyquest... Keep "fc + fm + dev" less then 40 percent or so of Fs or the
sidebands will start to fold over around Fs/2, and the signal is no longer
faithfully represented. Plots: with the "zoom on" enabled, you can take your
mouse into a figure and zoom in...
Regards,
bill schintler
fc4;
Fs*fc;
t=(0:Fs)'/Fs;
L0*fc;
fm0;
ph=0;
devH0;
NFFT=2^nextpow2(L);
f=Fs*linspace(0,1,NFFT);
%%---- FM ----%%
x=1*cos((2*pi*fc*t)-(dev/fm)*1*sin(2*pi*fm*t+(ph*pi/180)));
[num den]=butter(1,(fc+30+30)/(Fs/2)); % -3db @ 10Hz
x1 = filtfilt(num,den,x); % leave DC
%%---- FM DEMOD ----%%
yq = hilbert(x1).*exp(-j*2*fc*t);
y = (1/(2*pi*dev))*[zeros(1,size(yq,2)); diff(unwrap(angle(yq)))*Fs];
[num den]=butter(1,25/(Fs/2),'high');
y = filtfilt(num,den,y);
[num den]=butter(1,30/(Fs/2));
y = filtfilt(num,den,y);
figure;
subplot(211)
plot(t,y)
grid;
zoom on;
%window
Window = hanning(length(y));
subplot(212)
plot(f,2*abs(fft(Window.*y,NFFT)/L))
grid;
zoom on;
%axis([0 200 0 4])