基于matlab仿真的数字调制与解调设计

随着通信系统复杂性的增加，传统的手工分析与电路板试验等分析设计方法已经不能适应发展的需要，通信系统计算机模拟仿真技术日益显示出其巨大的优越性。计算机仿真是根据被研究的真实系统的模型，利用计算机进行实验研究的一种方法.它具有利用模型进行仿真的一系列优点，如费用低，易于进行真实系统难于实现的各种试验，以及易于实现完全相同条件下的重复试验等。Matlab仿真软件就是分析通信系统常用的工具之一。

Matlab是一种交互式的、以矩阵为基础的软件开发环境，它用于科学和工程的计算与可视化。Matlab的编程功能简单，并且很容易扩展和创造新的命令与函数。应用Matlab可方便地解决复杂数值计算问题。Matlab具有强大的Simulink动态仿真环境，可以实现可视化建模和多工作环境间文件互用和数据交换。Simulink支持连续、离散及两者混合的线性和非线性系统，也支持多种采样速率的多速率系统，Simulink为用户提供了用方框图进行建模的图形接口，它与传统的仿真软件包用差分方程和微分方程建模相比，更直观、方便和灵活。用户可以在Matlab和Simulink两种环境下对自己的模型进行仿真、分析和修改。用于实现通信仿真的通信工具包（Communication toolbox，也叫Commlib，通信工具箱）是Matlab语言中的一个科学性工具包，提供通信领域中计算、研究模拟发展、系统设计和分析的功能，可以在Matlab环境下独立使用，也可以配合Simulink使用。另外，Matlab的图形界面功能GUI（Graphical User Interface）能为仿真系统生成一个人机交互界面，便于仿真系统的操作。因此，Matlab在通信系统仿真中得到了广泛应用，本文也选用该工具对数字调制系统进行仿真。

1.2　数字通信的发展现状和趋势

进入20世纪以来，随着晶体管、集成电路的出现与普及、无线通信迅速发展。特别是在20世纪后半叶，随着人造地球卫星的发射，大规模集成电路、电子计算机和光导纤维等现代技术成果的问世，通信技术在以下几个不同方向都取得了巨大的成功。

（1）微波中继通信使长距离、大容量的通信成为了现实。

（2）移动通信和卫星通信的出现，使人们随时随地可通信的愿望可以实现。

（3）光导纤维的出现更是将通信容量提高到了以前无法想象的地步。

（4）电子计算机的出现将通信技术推上了更高的层次，借助现代电信网和计算机的融合，人们将世界变成了地球村。

（5）微电子技术的发展，使通信终端的体积越来越小，成本越来越低，范围越来越广。例如，2003年我国的移动电话用户首次超过了固定电话用户。根据国家信息产业部的统计数据，到2005年底移动电话用户近4亿。

随着现代电子技术的发展，通信技术正向着数字化、网络化、智能化和宽带化的方向发展。随着科学技术的进步，人们对通信的要求越来越高，各种技术会不断地应用于通信领域，各种新的通信业务将不断地被开发出来。到那时人们的生活将越来越离不开通信。

1.3　研究目的与意义

数字调制是指用数字基带信号对载波的某些参量进行控制，使载波的这些参量随基带信号的变化而变化。根据控制的载波参量的不同，数字调制有调幅、调相和调频三种基本形式，并可以派生出多种其他形式。由于传输失真、传输损耗以及保证带内特性的原因，基带信号不适合在各种信道上进行长距离传输。为了进行长途传输，必须对数字信号进行载波调制，将信号频谱搬移到高频处才能在信道中传输。因此，大部分现代通信系统都使用数字调制技术。另外，由于数字通信具有建网灵活，容易采用数字差错控制技术和数字加密，便于集成化，并能够进入综合业务数字网（ISDN网），所以通信系统都有由模拟方式向数字方式过渡的趋势。因此，对数字通信系统的分析与研究越来越重要，数字调制作为数字通信系统的重要部分之一，对它的研究也是有必要的。通过对调制系统的仿真，我们可以更加直观的了解数字调制系统的性能及影响性能的因素，从而便于改进系统，获得更佳的传输性能。

2PSK调制解调程序

clear;
t0=0.15;
ts=0.001;
fc=200;
snr=10;
fs=1/ts;
df=0.2;
t=[ts:ts:t0];
snr_lin=10^(snr/10);
y=[];
nm=length(t);
 
%生成2PSK
for i=1:10
   x=rand;
      for j=1:15
      l=(i-1)*15+j;
          c(l)=cos(2*pi*fc*t(l));
          c1(l)=cos(2*pi*fc*t(l)+pi);
      if x>=0.5
         mi(l)=1;
         y(l)=c(l);
      elseif x<0.5
         mi(l)=0;
         y(l)=c1(l);
      end
     end
end
%加噪声
signal_power=spower(y(1:length(t)));       % power in modulated signal
noise_power=signal_power/snr_lin;          % compute noise power
noise_std=sqrt(noise_power);              % compute noise standard deviation
noise=noise_std*randn(1,length(y));  % generate noise
r=y+noise;                              % add noise to the modulated signal
%时域频域转换
m=mi;
[M,m,df1]=fftseq(m,ts,df);
M=M*ts;
f=[0:df1:df1*(length(m)-1)]-fs/2;
[Y,y,df1]=fftseq(y,ts,df);
Y=Y*ts;
[C,c,df1]=fftseq(c,ts,df);
C=C*ts;
[R,r,df1]=fftseq(r,ts,df);
R=R*ts;
[NOISE,noise,df1]=fftseq(noise,ts,df);
NOISE=NOISE*ts;
%解调
f_cutoff=70;%70Hz低通滤波器
n_cutoff=floor(70/df1);
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);
yy=r.*c;
[YY,yy,df1]=fftseq(yy,ts,df);
YY=YY*ts;
DEM=H.*YY;%滤波
dem=real(ifft(DEM))/fs;
dem=dem(1:length(t));
yout=dem;
for i=(1:length(t));% 判决，得到解调结果          
    if yout(1,i)>0;    
        yout(1,i)=1;
    else
        yout(1,i)=0;
    end;
end;
 [YOUT,yout,df1]=fftseq(yout,ts,df);
YOUT=YOUT*ts;
pause;
%仿真
pause;
figure(1);
subplot(2,2,1);
plot(t,mi(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('源信号波形');
subplot(2,2,2);
plot(f,abs(fftshift(M)));xlabel('Frequency'); title('源信号频谱');
subplot(2,2,3);
plot(t,c(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('载波波形');
subplot(2,2,4);
plot(f,abs(fftshift(C)));xlabel('Frequency');title('载波频谱');
pause;
figure(2)
subplot(3,2,1);
plot(t,noise(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('噪声波形');
subplot(3,2,2);
plot(f,abs(fftshift(NOISE)));xlabel('Frequency');title('噪声频谱');
subplot(3,2,3);
plot(t,y(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('未加噪声调制波形');
subplot(3,2,4);
plot(f,abs(fftshift(Y)));xlabel('Frequency'); title('未加噪声调制频谱');
subplot(3,2,5);
plot(t,r(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('加噪声调制波形');
subplot(3,2,6);
plot(f,abs(fftshift(R)));xlabel('Frequency'); title('加噪声调制频谱');
pause;
figure(3);
subplot(2,2,1);
plot(t,yy(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('相干解调后波形');
subplot(2,2,2);
plot(f,abs(fftshift(YY)));xlabel('Frequency'); title('相干解调后频谱');
subplot(2,2,3);
plot(t,dem(1:length(t)));grid;xlabel('Time');title('低通后波形');
subplot(2,2,4);
plot(f,abs(fftshift(DEM)));xlabel('Frequency'); title('低通后频谱');
pause;
figure(4)
subplot(2,2,1);
plot(t,yout(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('抽样判决后波形');
subplot(2,2,2);
plot(f,abs(fftshift(YOUT)));xlabel('Frequency'); title('抽样判决后频谱');
subplot(2,2,3);
plot(t,mi(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('源信号波形');
subplot(2,2,4);
plot(f,abs(fftshift(M)));xlabel('Frequency'); title('源信号频谱');
 
2DPSK调制解调程序
clear;
t0=0.15;
ts=0.001;
fc=200;
snr=10;
fs=1/ts;
df=0.2;
t=[ts:ts:t0];
snr_lin=10^(snr/10);
y=[];
nm=length(t);
%生成2DPSK
bn=[];bn(1)=1;%bn是差分码，设bn的第一个符号为1
for i=1:10
   x=rand;
   if x>=0.5
      m(i)=1;
    else
      m(i)=0;   
   end
   bn(i+1)=xor(m(i),bn(i));
   for j=1:15
        l=(i-1)*15+j;
        c(l)=cos(2*pi*fc*t(l));
        c1(l)=cos(2*pi*fc*t(l)+pi);
      if m(i)==1;
         mi(l)=1;
      else
         mi(l)=0;
      end
      if bn(i+1)==1;
         bn1(l)=1;
         y(l)=c(l);
      else
         bn1(l)=0;
         y(l)=c1(l);
      end
   end
end
%加噪声
signal_power=spower(y(1:length(t)));       % power in modulated signal
noise_power=signal_power/snr_lin;          % compute noise power
noise_std=sqrt(noise_power);              % compute noise standard deviation
noise=noise_std*randn(1,length(y));  % generate noise
r=y+noise;                              % add noise to the modulated signal
%时域频域转换
m=mi;
[M,m,df1]=fftseq(m,ts,df);
M=M*ts;
[BN1,bn1,df1]=fftseq(bn1,ts,df);
BN1=BN1*ts;
f=[0:df1:df1*(length(m)-1)]-fs/2;
[Y,y,df1]=fftseq(y,ts,df);
Y=Y*ts;
[C,c,df1]=fftseq(c,ts,df);
C=C*ts;
[R,r,df1]=fftseq(r,ts,df);
R=R*ts;
[NOISE,noise,df1]=fftseq(noise,ts,df);
NOISE=NOISE*ts;
%解调
f_cutoff=70;
n_cutoff=floor(70/df1);
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);
yy=r.*c;
[YY,yy,df1]=fftseq(yy,ts,df);
YY=YY*ts;
DEM=H.*YY;%滤波
dem=real(ifft(DEM))/fs;
dem=dem(1:length(t));
ytemp=dem;
y1=ytemp;
for i=(1:length(t));% 判决，得到解调结果          
    if y1(1,i)>0;    
       y1(1,i)=1;
    else
        y1(1,i)=0;
    end;
end;
for i=(1:10);%码反变换
    k=y1((i-1)*15+1:i*15);
    if mean(k)>0.5;
       y2(i)=1;    
     else
       y2(i)=0;
    end  
end
y3(1)=xor(y2(1),1);
for i=(2:10);
    y3(i)=xor(y2(i-1),y2(i));
end
for i=(1:10);
    for j=(1:15);
        l=(i-1)*15+j;
        if y3(i)==1;
           yout(l)=1;
        else
           yout(l)=0;
        end
    end 
end 
[Y1,y1,df1]=fftseq(y1,ts,df);
Y1=Y1*ts;
[YOUT,yout,df1]=fftseq(yout,ts,df);
YOUT=YOUT*ts;
pause;
%仿真
pause;
figure(1);
subplot(3,2,1);
plot(t,mi(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('源信号波形');
subplot(3,2,2);
plot(f,abs(fftshift(M)));xlabel('Frequency'); title('源信号频谱');
subplot(3,2,3);
plot(t,bn1(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('码变换后波形');
subplot(3,2,4);
plot(f,abs(fftshift(BN1)));xlabel('Frequency');title('码变换后频谱');
subplot(3,2,5);
plot(t,c(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('载波波形');
subplot(3,2,6);
plot(f,abs(fftshift(C)));xlabel('Frequency');title('载波频谱');
pause;
figure(2);
subplot(3,2,1);
plot(t,noise(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('噪声波形');
subplot(3,2,2);
plot(f,abs(fftshift(NOISE)));xlabel('Frequency');title('噪声频谱');
subplot(3,2,3);
plot(t,y(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('未加噪声调制波形');
subplot(3,2,4);
plot(f,abs(fftshift(Y)));xlabel('Frequency'); title('未加噪声调制频谱');
subplot(3,2,5);
plot(t,r(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('加噪声调制波形');
subplot(3,2,6);
plot(f,abs(fftshift(R)));xlabel('Frequency'); title('加噪声调制频谱');
pause;
figure(3);
subplot(3,2,1);
plot(t,yy(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('相干解调后波形');
subplot(3,2,2);
plot(f,abs(fftshift(YY)));xlabel('Frequency'); title('相干解调后频谱');
subplot(3,2,3);
plot(t,dem(1:length(t)));grid;xlabel('Time');title('低通后波形');
subplot(3,2,4);
plot(f,abs(fftshift(DEM)));xlabel('Frequency'); title('低通后频谱');
subplot(3,2,5);
plot(t,y1(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('抽样判决后波形');
subplot(3,2,6);
plot(f,abs(fftshift(Y1)));xlabel('Frequency'); title('抽样判决后频谱');
pause;
figure(4);
subplot(2,2,1);
plot(t,yout(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('码反变换后波形');
subplot(2,2,2);
plot(f,abs(fftshift(YOUT)));xlabel('Frequency'); title('码反变换后频谱');
subplot(2,2,3);
plot(t,mi(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('源信号波形');
subplot(2,2,4);
plot(f,abs(fftshift(M)));xlabel('Frequency'); title('源信号频谱');
 
16QAM调制解调程序
clear;
t0=0.15;
ts=0.001;
fc=200;
snr=10;
fs=1/ts;
df=0.2;
t=[ts:ts:t0];
snr_lin=10^(snr/10);
nm=length(t);
for i=1:10
   x=rand;
   for j=1:15
       l=(i-1)*15+j;
       c(l)=cos(2*pi*fc*t(l));
       s(l)=sin(2*pi*fc*t(l));
      if x>=0.7
         mi(l)=3;
      elseif x<=0.3
         mi(l)=1;
      elseif x<=0.5
          mi(l)=-1;
      else
          mi(l)=-3;
      end
   end
end
for i=1:10
   x=rand;
   for j=1:15
       l=(i-1)*15+j;
       c(l)=cos(2*pi*fc*t(l));
       s(l)=sin(2*pi*fc*t(l));
      if x>=0.8
         mq(l)=3;
      elseif x<=0.3
         mq(l)=1;
      elseif x<=0.6
          mq(l)=-1;
      else
          mq(l)=-3;
      end
  end
end
nm=length(t);
mi(nm)=mi(nm-1);
mq(nm)=mq(nm-1);
ask1=mi.*c;
ask2=mq.*s;
y=mi.*c+mq.*s;
%加噪声
signal_power=spower(y(1:length(t)));       % power in modulated signal
noise_power=signal_power/snr_lin;          % compute noise power
noise_std=sqrt(noise_power);              % compute noise standard deviation
noise=noise_std*randn(1,length(y));  % generate noise
r=y+noise;                              % add noise to the modulated signal
yi=r.*c;
yq=r.*s;
%%时域频域转换;
mk=mi;
ml=mq;
[M1,mk,df1]=fftseq(mk,ts,df);
M1=M1/fs;
[M2,ml,df1]=fftseq(ml,ts,df);
M2=M2/fs;
[C,c,df1]=fftseq(c,ts,df);
C=C/fs;
[S,s,df1]=fftseq(s,ts,df);
S=S/fs;
[ASK1,ask1,df1]=fftseq(ask1,ts,df);
ASK1=ASK1/fs;
[ASK2,ask2,df1]=fftseq(ask2,ts,df);
ASK2=ASK2/fs;
f=[0:df1:df1*(length(mk)-1)]-fs/2;
[NOISE,noise,df1]=fftseq(noise,ts,df);
NOISE=NOISE*ts;
u=y;
[U,u,df1]=fftseq(u,ts,df);
U=U/fs;
[R,r,df1]=fftseq(r,ts,df);
R=R/fs;
[Yi,yi,df1]=fftseq(yi,ts,df);
Yi=Yi/fs;
[Yq,yq,df1]=fftseq(yq,ts,df);
Yq=Yq/fs;
%滤波
f_cutoff=70;
n_cutoff=floor(70/df1);
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);
DEM1=H.*Yi;
dem1=real(ifft(DEM1))/fs;
DEM2=H.*Yq;
dem2=real(ifft(DEM2))/fs;
dem1=dem1(1:length(t));
dem2=dem2(1:length(t));
%抽样判决
for i=(1:10);
    k=dem1((i-1)*15+1:i*15);
    x1(i)=mean(k);
    for j=(1:15)
        l=(i-1)*15+j;
      if x1(i)>=0
          x1(i)=max(k)
          if x1(i)>3*10^-6
             y1(l)=3
          else
             y1(l)=1
          end
      end
      if x1(i)<0
          x1(i)=min(k)
          if x1(i)<-3*10^-6
             y1(l)=-3
          else
              y1(l)=-1
          end
      end
    end 
end
for i=(1:10);
    k=dem2((i-1)*15+1:i*15);
    x2(i)=mean(k);
    for j=(1:15)
        l=(i-1)*15+j;
      if x2(i)>=0
          x2(i)=max(k)
          if x2(i)>3*10^-6
             y2(l)=3
          else
             y2(l)=1
          end
      end
      if x2(i)<0
          x2(i)=min(k)
          if x2(i)<-3*10^-6
             y2(l)=-3
          else
             y2(l)=-1
          end
      end
    end 
end
 [Y1,y1,df1]=fftseq(y1,ts,df);
Y1=Y1/fs;
[Y2,y2,df1]=fftseq(y2,ts,df);
Y2=Y2/fs;
%仿真
pause;
figure(1)
subplot(2,2,1)
plot(t,mi(1:length(t)));grid;axis([0 0.15 -3.5 3.5]);xlabel('Time');title('源信号1波形');
subplot(2,2,2)
plot(f,abs(fftshift(M1)));xlabel('Frequncy');title('源信号1频谱');
subplot(2,2,3)
plot(t,mq(1:length(t)));grid;axis([0 0.15 -3.5 3.5]);xlabel('Time');title('源信号2波形');
subplot(2,2,4)
plot(f,abs(fftshift(M2)));xlabel('Frequncy');title('源信号2频谱');
pause;
figure(2)
subplot(2,2,1)
plot(t,c(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('载波1波形');
subplot(2,2,2)
plot(f,abs(fftshift(C)));xlabel('Frequncy');title('载波1频谱');
subplot(2,2,3)
plot(t,s(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('载波2波形');
subplot(2,2,4)
plot(f,abs(fftshift(S)));xlabel('Frequncy');title('载波2频谱');
pause;
figure(3)
subplot(2,2,1)
plot(t,ask1(1:length(t)));grid;axis([0 0.15 -5 5]);xlabel('Time');title('4ASK信号1波形');
subplot(2,2,2)
plot(f,abs(fftshift(ASK1)));xlabel('Frequncy');title('4ASK信号1频谱');
subplot(2,2,3)
plot(t,ask2(1:length(t)));grid;axis([0 0.15 -5 5]);xlabel('Time');title('4ASK信号2波形');
subplot(2,2,4)
plot(f,abs(fftshift(ASK1)));xlabel('Frequncy');title('4ASK信号2频谱');
pause;
figure(4)
subplot(3,2,1)
plot(t,noise(1:length(t)));grid;axis([0 0.15 -2 2]);xlabel('Time');title('噪声波形');
subplot(3,2,2)
plot(f,abs(fftshift(NOISE)));xlabel('Frequency');title('噪声频谱');
subplot(3,2,3)
plot(t,u(1:length(t)));grid;axis([0 0.15 -5 5]);xlabel('Time');title('未加噪声调制16QAM波形');
subplot(3,2,4)
plot(f,abs(fftshift(U)));xlabel('Frequncy');title('未加噪声调制16QAM频谱');
subplot(3,2,5)
plot(t,r(1:length(t)));grid;axis([0 0.15 -5 5]);xlabel('Time');title('加噪声调制16QAM波形');
subplot(3,2,6)
plot(f,abs(fftshift(R)));xlabel('Frequncy');title('加噪声调制16QAM频谱');
pause;
figure(5)
subplot(2,2,1)
plot(t,yi(1:length(t)));grid;xlabel('Time');title('相干解调后信号1波形');
subplot(2,2,2)
plot(f,abs(fftshift(Yi)));xlabel('Frequncy');title('相干解调后信号1频谱');
subplot(2,2,3)
plot(t,dem1(1:length(t)));grid;xlabel('Time');title('低通后信号1波形');
subplot(2,2,4)
plot(f,abs(fftshift(DEM1)));xlabel('Frequency');title('低通后信号1频谱');
pause;
figure(6)
subplot(2,2,1)
plot(t,yq(1:length(t)));grid;xlabel('Time');title('相干解调后信号2波形');
subplot(2,2,2)
plot(f,abs(fftshift(Yq)));xlabel('Frequency');title('相干解调后信号2频谱');
subplot(2,2,3)
plot(t,dem2(1:length(t)));grid;xlabel('Time');title('低通后信号2波形');
subplot(2,2,4)
plot(f,abs(fftshift(DEM2)));xlabel('Frequency');title('低通后信号2频谱');
pause;
figure(7)
subplot(2,2,1)
plot(t,y1(1:length(t)));grid;axis([0 0.15 -3.5 3.5]);xlabel('Time');title('抽样判决后信号1波形');
subplot(2,2,2)
plot(f,abs(fftshift(Y1)));xlabel('Frequency');title('抽样判决后信号1频谱');
subplot(2,2,3)
plot(t,mi(1:length(t)));grid;axis([0 0.15 -3.5 3.5]);xlabel('Time');title('源信号1波形');
subplot(2,2,4)
plot(f,abs(fftshift(M1)));xlabel('Frequncy');title('源信号1频谱');
pause;
figure(8)
subplot(2,2,1)
plot(t,y2(1:length(t)));grid;axis([0 0.15 -3.5 3.5]);xlabel('Time');title('抽样判决后信号2波形');
subplot(2,2,2)
plot(f,abs(fftshift(Y2)));xlabel('Frequency');title('抽样判决后信号2频谱')
subplot(2,2,3)
plot(t,mq(1:length(t)));grid;axis([0 0.15 -3.5 3.5]);xlabel('Time');title('源信号2波形');
subplot(2,2,4)
plot(f,abs(fftshift(M2)));xlabel('Frequncy');title('源信号2频谱');
 
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;
 
 
spower.m程序
function p=spower(x)
%       p=spower(x)
%SPOWER     returns the power in signal x
p=(norm(x)^2)/length(x);