fc = 150;                     %載波頻率
fs = 400;                     %采樣頻率
t = (0 :0.001:0.15);             %時間區域
x = sin(2*pi*30*t);             %輸入信號
y = modulate(x,fc,fs,'FM');       %調制信號
plot(t,y);                     %FM信號時域圖
xlabel('t(s)');ylabel('y');
axis([0 0.15 -1 1]);
b=fft(y,1024);                 %對y進行傅利葉變換
f=(0:length(b)-1)*fs/length(b) -fs/2;
figure
plot(f,abs(b));                 %FM信號頻譜圖
xlabel('Frequence(Hz)');         %頻率分量
ylabel('Power Spectrum(y)');     %頻率譜能量