振幅调制和解调的原理及MATLAB编程实现 已给出一个编码,但有错误,求修改!

振幅调制和解调的原理及MATLAB编程实现 已给出一个编码,但有错误,求修改!
% MATLAB script for Illustrative Problem 3.5.
% Demonstration script for DSB-AM demodulation.The message signal
% is +1 for 0 < t < t0/3,-2 for t0/3 < t < 2t0/3,and zero otherwise.
echo on
t0=.15; % 信号的持续时间
ts=1/1500; % 采样间隔
fc=250; % 载波频率
fs=1/ts; % 采样频率
t=[0:ts:t0]; % 时间向量
df=0.3; % 所需的频率分辨率
% 信息信号
m=[ones(1,t0/(3*ts)),-2*ones(1,t0/(3*ts)),zeros(1,t0/(3*ts)+1)];
c=cos(2*pi*fc.*t); %
u= ones(1,length(c) )
y= ones(1,length(c) )
[M,m,df1]=fftseq(m,ts,df); %
M=M/fs; % 缩放
[U,u,df1]=fftseq(u,ts,df); %
U=U/fs; % 缩放
[Y,y,df1]=fftseq(y,ts,df); %
Y=Y/fs; % 缩放
f_cutoff=150; % cutoff freq.of the filter
n_cutoff=floor(150/df1); % Design the filter.
f=[0:df1:df1*(length(y)-1)]-fs/2;
H=zeros(size(f));
H(1:n_cutoff)=2*ones(1,n_cutoff);
H(length(f)-n_cutoff+1:length(f))=2*ones(1,n_cutoff);
DEM=H.*Y; % spectrum of the filter output
dem=real(ifft(DEM))*fs; % filter output
pause % Press a key to see the effect of mixing.
clf
subplot(3,1,1)
plot(f,fftshift(abs(M)))
title('信息信号的频谱')
xlabel('频率')
subplot(3,1,2)
plot(f,fftshift(abs(U)))
title('调制信号的频谱')
xlabel('频率')
subplot(3,1,3)
plot(f,fftshift(abs(Y)))
title('混频器的输出频谱')
xlabel('频率')
pause % Press a key to see the effect of filtering on the mixer output.
clf
subplot(3,1,1)
plot(f,fftshift(abs(Y)))
title('混频器的输出频谱')
xlabel('Frequency')
subplot(3,1,2)
plot(f,fftshift(abs(H)))
title('低通滤波特性')
xlabel('Frequency')
subplot(3,1,3)
plot(f,fftshift(abs(DEM)))
title('解调器的输出频谱')
xlabel('Frequency')
pause % Press a key to compare the spectra of the message and the received signal.
clf
subplot(2,1,1)
plot(f,fftshift(abs(M)))
title('信息信号的频谱')
xlabel('Frequency')
subplot(2,1,2)
plot(f,fftshift(abs(DEM)))
title('解调器的输出频谱')
xlabel('Frequency')
pause % Press a key to see the message and the demodulator output signals.
subplot(2,1,1)
plot(t,m(1:length(t)))
title('信息信号')
xlabel('Time')
subplot(2,1,2)
plot(t,dem(1:length(t)))
title('该解调器输出')
xlabel('Time')

补上所缺函数fftseq.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;