echo ont0=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