oversampling_echo_gen.m

天气雷达回波产生 这个代码是自己编的 我用过没什么问题

function [oversampling_echo]=oversampling_echo_gen(M,L,h,p,fs,v,dv,lamga,Ps,Pn)
%M: 时间采样点数 
%L：过采样率
%h: 滤波器冲击响应
%p: 脉冲包络系数
%fs:重复频率
%v：信号多普勒速度
%dv：信号多普勒速度谱宽
%lamga：雷达波长
%Ps：信号功率，单位dbm
%Pn：噪声功率，单位dbm
%oversampling_echo  输出  M x L

% close all;
% clc;


 M=128;
% %range 
 L=10;
 N=3*L;
% %oversampling rate 

%transmitter pulse shape
 for i=1:L
    p(i)=1;
end
%receiver filter response
%filter response length

h=[0.077630037708532,0.189678324294802,0.306045767992654,0.40049650976237,0.450257096821268];
F=length(h);
N=N+2*F-1+L/2;
%refrectivity gradient 
for i=1:L
    g(i)=1;
end
%construct MxL sample time data
%信号表达式：
%V(1,q)=[leijia(s(l+i).*g(i,n).*p(L-1-i))]*h(1)
%l=0,1...,L+F-2;n=0,1,..M-1;

%********************************************************************
%************* 构造range 相关*****************************************
%********************************************************************
% 一个径向M个radia,一个radia共L个数据
%产生M x N的zero mean,variance one random
SS=randn(M,N);
for x=1:M
%随机相位
a=2*pi*randn(1,N);
b=SS(1,:);
%目标回波
S=b.*(sin(a)+j*cos(a));
for i=1:M
SC(i,:)=SS(i,:).*(sin(a)+j*cos(a));
end
S1=zeros(1,N-L);
for l=1:N-L
    for i=1:L
        S1(l)=S(l+i-1)*p(L-i+1)*g(i)+S1(l);  % 累加
    end
end
VC=conv(S1,h);
%截取
VQ=VC(F:length(VC)-F+1);
%产生不同 radia的数据
%%%屏蔽h的作用
VP(x,:)=VQ;
%VP(x,:)=S1;
end
 figure
 subplot(311);plot(real(VP(:,1)));title('均值为0，方差为1的高斯随机信号');
xx=VP(1,:);
xxx=xx(1:L);
 yy=S(1,:);
 yyy=yy(1:L);
 subplot(312);plot(real(xcorr(xxx)));title('粒子叠加高斯随机信号自相关函数');
subplot(313);plot(real(xcorr(yyy)));title('独立高斯随机信号自相关函数');
%重复频率
%****************************************************************
%*************** 构造sample 相关*********************************
%****************************************************************
%%%设计参数
 fs=1270;
 lamga=0.1;
G=sum(conv(p,h).^2);
G=8;
 v=10;
f=2*v/lamga;
 dv=2;
df=2*dv/lamga;
SNR=-20;
cc=0.00000000;
Pnn=fs*cc;
Pr=Pnn*10^(SNR/10);
VA=lamga*fs/2;
FA=fs/2;
%**高斯谱模型
k=-M:2*M-1;
VK=-VA+(2*k/M)*VA;
 FK=-FA+(2*k/M)*FA;
FK=-2*VK/lamga;
%如何归一化?
NPs=128;
VVK=(Pr*1/(G*sqrt(2*pi)*dv))*exp(-(VK-v).^2/(2*dv^2))+cc;
VFK=(1/(sqrt(2*pi)*df))*exp(-(FK-f).^2/(2*df^2))+cc;
%  figure,
 plot(k,VFK);title('频率谱（单位：Hz)');
%%截断为Nyquist频率(-VA,VA)(-FS/2,FS/2)
SIG_V=sqrt(VVK(M:2*M-1));
%%%对功率谱进行规一化，方法：除以平均功率
SIG_F=sqrt(VFK(M:2*M-1)/mean(VFK(M:2*M-1)));
nn_vk=FK(M:2*M-1);
 figure,
plot(nn_vk,20*log(SIG_F));title('频率谱（单位：Hz)');

%***************************************************************
%************** sample 相关， range 不相关 ********************
%*************************************************************** 
[ROW_SC,COL_SC]=size(SC);
for i=1:COL_SC
    aa=(fft(SC(:,i)));
    %%%对功率谱进行规一化，方法：除以平均功率
    bb=aa/mean(abs(aa).^2);
    SQ(:,i)=(ifft(aa.*fftshift(SIG_F'))); 
end
nn=1:M;
    figure,
    subplot(211);plot(nn,real(SQ(:,3)),'r',nn,imag(SQ(:,3)),'b');
    title('加噪声前高斯谱模型天气回波的时域波形');
    subplot(212);plot(nn_vk,10*log10(abs(fftshift(fft(SQ(:,3)))))','b');
    title('对应频谱分布（单位：Hz）');
%***************************************************************
%************** sample and range 相关合二为1 ********************
%*************************************************************** 

%%% 注意：做ifft时，要先fftshift
[ROW_VP,COL_VP]=size(VP);
for i=1:COL_VP
    aa=(fft(VP(:,i)));
    bb=0.001.*aa/mean(abs(aa).^2);
   xx=aa.*fftshift(SIG_F');
  figure,subplot(311);plot(10*log10(abs((real(bb)))));subplot(312);plot(10*log10(abs((SIG_F'))));subplot(313);plot(10*log10(abs(bb.*SIG_F')));
    XXQ(:,i)=((ifft(aa.*fftshift(SIG_F')))); 
    % 除平均功率归一化
    XQ(:,i)=XXQ(:,i)/sqrt(mean(abs(XXQ(:,i)).^2));  
end
    nn=1:M;
   figure,
     subplot(211);plot(nn,real(XQ(:,3)),'r',nn,imag(XQ(:,3)),'b');
     title('加噪声前高斯谱模型天气回波的时域波形');
     subplot(212);plot(nn_vk,10*log10(abs((fftshift(fft(XQ(:,1))))))','b');
     title('对应频谱分布（单位：Hz）');
% %**************************************************************
%*************** 时域 add noise 加入噪声****************************
%*************************************************************

%  -130dbm=1.00E-4mVpp,-100dbm=3.16E-3mVpp,-50dbm=1.00EmVpp
XTQ=XQ;
[ROW_VP,COL_VP]=size(XTQ);
Ps=15;%dBm
%dBm=10*log(((Vpp*10^-3)*0.707)^2)*1000/50);
 Pn=-80;%-50dBm = 1.00EmVpp
NEF=sqrt(50*(10^(Pn/10)/1000))*1.414;
 Ps=-95;%-30dBm = 1.00EmVpp
Pr=sqrt(50*(10^(Ps/10)/1000))*1.414;;
    XNQ=Pr*XTQ+ NEF*(randn(ROW_VP,COL_VP)+sqrt(-1)*randn(ROW_VP,COL_VP));
    oversampling_echo=XNQ(:,(1:L));
  oversampling_echo=XNQ;
  nn_vk=-fs/2:fs/M:fs/2-fs/M;
    figure(2);
    subplot(311);plot(nn,real(XNQ(:,3)),'r',nn,imag(XNQ(:,3)),'b');
     title('加噪声后的高斯谱模型天气回波的时域波形（单位：mv）');
   subplot(312);plot(nn_vk,20*log10((abs(fftshift(fft(XNQ(:,3)))))));title('频率谱（单位：Hz)');
%  yyy=XNQ(StartP:L+StartP-1);
    n_cor=1:2*L-1;
    subplot(313);plot(n_cor,real(xcorr(yyy)),'r',n_cor,imag(xcorr(yyy)),'b');title('过采样信号的自相关函数');