谱相关估计-频域平滑算法及matlab仿真

       对于频域平滑算法的估计式子

 = （1）

       设在实际接收的数据长度为Δt，采样周期为进行均匀抽样得到N个点，要求N为2的整数次幂的离散序列x(t,n),显然有。对离散序列x(t,n)做DFT，得到：

（2）

       上式中为采样间隔（亦称为频率分辨率）

       设频率平滑间隔为，M为参与平滑的样本数，则将（1）式离散化，用求和代替积分得到：

              (3)

       对于离散频率情况：设为谱频率，为循环频率，k和h分别为频率f和循环频率α基于采样间隔Fs的数字频率，则3式改写为：

（4）

算法流程图如下所示：

 

%用于求BPSK调制信号的循环谱  参考《基于频域平滑循环周期图法的
% 直接序列扩频信号的参数估计》
tic
clc
clear ;
clear all;
%%%%%%%%%%%%%%%仿真参数设置%%%%%%%%%%%%%%%%%
fd=1;                                      %数据速率为1k
Ts=1/fd;                                   %码元间隔Ts=0.001s
delt=2^7;                                  %时间段 单位是ms
fc=2;                                      %载波频率2k
fs=8;                                      %采样频率8k
dt=1/fs;                                   %采样时间间隔Ts=1/fs
num=fd*delt;                               %码元个数
t=[0:dt:delt-dt];                          %时间向量
%%%%%%%%%%%%%%产生数据%%%%%%%%%%%%%%%%%%%%%%     
data=sign(randn(1,num));
%%%%%%通过基带成型滤波器后的信号(矩形成型)%%%%%%%%%%%% 矩形信号成形就是线性插值？
 for rr=1:num
    I((rr-1)*fs+1:rr*fs)=data(rr);
end
%%%%%%%最终BPSK调制信号%%%%%%%%%%%%%%%%%%%%%
xx=I.*cos(1*2*pi*fc.*t);

%x=awgn(xx,-3);                             %信号加入高斯白噪声
%noise=wgn(1,1024,10);

s=fft(xx);                                 %求数据段的傅立叶变换
length(s)
F=fftshift(s);                             %将零频率搬至中心
Fc=conj(F);                                %取共轭
ff=fs.*(-length(xx)/2:length(xx)/2-1)/length(xx);      % 原信号的模拟频率f=fs*k/N
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=128;                                     %频率平滑点数
Fs=1;                                      %频率分辨率 fs/N 参考《基于循环平稳特征的频谱感知技术研究》（2-25）对Fs的定义
y=length(xx)-M/2;
%%%%%%计算偶数时间段%%%%%%
ress=zeros(length(xx)-M/2-M,2*(length(xx)-M/2)-M);
for i=M+1:length(xx)-M/2 % i为谱频率
  for k=-y:2:y-1 % k就是循环频率参考文献中公式3的α
      if and(and((i+k/2>=M+1),(i+k/2<=length(xx)-M/2)),and((i-k/2>=M+1),(i-k/2)<=length(xx)-M/2))
         length(F(i+k/2+(-M/2:M/2)*Fs));
         length(Fc(i-k/2+(-M/2:M/2)*Fs));
         ress(i-M,k+2+y)=(1/M)*sum((1/delt)*F(i+k/2+(-M/2:M/2)*Fs).*Fc(i-k/2+(-M/2:M/2)*Fs));
      end
  end
end
%%%%%%计算奇数时移段%%%%%%%
for i=M+1.5:length(xx)-M/2+0.5
  for k=-y+1:2:y
      if and(and((i+k/2>=M+1),(i+k/2<=length(xx)-M/2)),and((i-k/2>=M+1),(i-k/2)<=length(xx)-M/2))
         length(F(i+k/2+(-M/2:M/2)*Fs));
         length(Fc(i-k/2+(-M/2:M/2)*Fs));
         ress(i-M-0.5,k+y)=(1/M)*sum((1/delt)*F(i+k/2+(-M/2:M/2)*Fs).*Fc(i-k/2+(-M/2:M/2)*Fs));
      end
  end
end
%%%%%%%%%%%%%%计算f矩阵%%%%%%%%%%%%%%%%%%%%%
e=[M+1-length(xx)/2:length(xx)-M/2-length(xx)/2]';
f=repmat(e,1,2*(length(xx)-M/2)-M).*fs/length(xx);
%%%%%%%%%%%%%%%计算alfa矩阵%%%%%%%%%%%%%%%%%
ee=[-(length(xx)-M/2):length(xx)-M/2-M-1];
alfa=repmat(ee,length(xx)-M/2-M,1).*fs/length(xx);
%%%%%%%%%%%%%%%绘图%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
S=abs(ress)/max(max(abs(ress)));
mesh(alfa,f,S);
ylabel('Frequency (Hz)');xlabel('Cycle frequency (Hz)'); 
title('bpsk循环谱128dianpinghua三维图');