http://s10/middle/59a069d5n9eb77287ffb9&690
             figure 1. 调制信号和调频信号

http://s15/middle/59a069d5n9eb77b327b4e&690
        figure 2. 调制信号频谱和调频信号频谱

% 主程序 频率调制

% fm1.m
% Matlab demonstration script for frequency modulation.The message signal
% is +2 for 0<t<t0/3,-2 for t0/3<t<2t0/3,and zero otherwise.
echo on
t0=0.15; %signal duration
ts=0.0005; %sample interval
fc=200; %载波 frequency
kf=50; %modulation index
fs=1/ts; %sampling frequency
t=[0:ts:t0]; %time vector
df=0.25; %required frequency resolution
% message signal
m=[2*ones(1,t0/(3*ts)),-2*ones(1,t0/(3*ts)),zeros(1,t0/(3*ts)+1)];
int_m(1)=0;
for i=1:length(t)-1 %integral of m
    int_m(i+1)=int_m(i)+m(i)*ts;
    echo off;
end
echo on;
[M,m,df1]=fftseq(m,ts,df); %Fourier transform
M=M/fs; %scaling
f=[0:df1:df1*(length(m)-1)]-fs/2; %frequency vector
u=cos(2*pi*fc*t+2*pi*kf*int_m); %modulated signal
[U,u,df]=fftseq(u,ts,df); %Fourier transform
U=U/fs; %scaling

pause %Press any key to see plot of the message and the modulated signal
subplot(211)
plot(t,m(1:length(t)))
axis([0 0.15 -2.1 2.1])
xlabel('Time')
title('The message signal')
subplot(212)
plot(t,u(1:length(t)))
axis([0 0.15 -2.1 2.1])
xlabel('Time')
title('The modulated signal')

pause %Press any key to see plot of the message and the modulated signal
subplot(211)
plot(f,abs(fftshift(M)))
xlabel('Frequency')
title('Magnitude spectrum of the message signal')
subplot(212)
plot(f,abs(fftshift(U)))
xlabel('Frequency')
title('Magnitude spectrum of the modulated signal')

%------------------------------------------------------------------

% -----主程序调用的函数，单独编程，程序名是MATLAB自动命名为函数名----

% MATLAB函数fftseq.m，将时间序列m，采样间隔ts和要求的频率分辨率df作为输入，就得长度为2的幂的序列，
% 这个序列的FFT(M)和所要的频率分辨率。
function [M,m,df]=fftseq(m,ts,df)
% [M,m,df]=fftseq(m,ts,df)
% [M,m,df]=fftseq(m,ts)
% FFTSEQ Generates M,the FFT of the sequence m.
% The sequence is zero-padded to meet the required frequency resolution df.
% ts is the sampling interval.The output df is the final frequency resolution.
% Output m is the zero-padded version of input m,M is the FFT.
fs=1/ts;
if nargin == 2
    n1=0;
else
    n1=fs/df;
end
n2=length(m);
n=2^(max(nextpow2(n1),nextpow2(n2)));
M=fft(m,n);
m=[m,zeros(1,n-n2)];
df=fs/n;
end

% 参考文献：约翰.G.普罗克斯  刘树棠译；现代通信系统使用MATLAB[M].西安交通大学出社,2001年，1版

% 程序清单来源 《现代通信系统 使用MATLAB》p118,p23

%（约翰.G.普罗克斯  刘树棠译 西安交通大学出社） 2001年11月 第1版