tdd_to_fdd6.m

不同频率条件下移动终端智能天线信噪比性能分析。。。。。。。。

%%%%%%%%%上行波束赋形
clear,clc;

M=4;             %array_number;
fc=2e9;          %work_frequency;
c=3e8;           %rate of light
lamda=c/fc;      %work_wavelength;
d=lamda/2;       %array_distance;

SNR=0;
JNR1=30;
JNR2=30;
JNR3=30;

ang0=15;
ang1=25;
ang2=0;
ang3=-15;

K=100;           %snap_number;
fs=1.8e10;       %sample_frequency;

a0=zeros(M,1);        % a0 期望信号方向矢量
a1=zeros(M,1);        % a1 干扰信号1的方向矢量
a2=zeros(M,1);        % a2 干扰信号2的方向矢量
a3=zeros(M,1);        % a3 干扰信号3的方向矢量
j=sqrt(-1);
for i=1:M
    a0(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang0*pi/180));
    a1(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang1*pi/180));
    a2(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang2*pi/180));
    a3(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang3*pi/180));
end
% A=[a0';a1';a2';a3']';
% A1=[a1';a2';a3']';
A=[a0,a1,a2,a3];
A1=[a1,a2,a3];

t=(1:1:K)/fs; 

s0=cos(2*pi*fc*t);        

AM_JAM1=1;
PM_JAM1=0;
AM_JAM2=1;
PM_JAM2=0;
AM_JAM3=1;
PM_JAM3=0;


if AM_JAM1==1
    s01=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM1==1
    s01=sin(2*pi*fc*t+randn(1,K));
end

if AM_JAM2==1
    s02=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM2==1
    s02=sin(2*pi*fc*t+randn(1,K));
end

if AM_JAM3==1
    s03=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM3==1
    s03=sin(2*pi*fc*t+randn(1,K));
end


Ps0=s0*s0'/K;                       % 计算期望信号的功率
Ps1=s01*s01'/K;                     % 计算干扰信号1的功率
Ps2=s02*s02'/K;                     % 计算干扰信号2的功率
Ps3=s03*s03'/K;                     % 计算干扰信号2的功率

Ps=10.^(SNR/10);
Pj1=10.^(JNR1/10);
Pj2=10.^(JNR2/10);
Pj3=10.^(JNR3/10);

s=s0/Ps0^0.5*Ps^0.5;             % 具有输入信噪比的期望信号
s1=s01/Ps1^0.5*Pj1^0.5;          % 具有干噪比1的干扰信号1
s2=s02/Ps2^0.5*Pj2^0.5;          % 具有干噪比2的干扰信号2
s3=s03/Ps3^0.5*Pj3^0.5;          % 具有干噪比2的干扰信号2

S=[s;s1;s2;s3];
S1=[s1;s2;s3];
numJam=size(S1,1);
Noise=randn(M,K);
X=A*S+Noise;
X1=A1*S1+Noise;

R=X*X'/K;
Rjn=X1*X1'/K;


h0=a0/(a0'*a0)^0.5;
B0=zeros(M-1,M);
for i=1:M-1
    for k=i:i+1
        B0(i,k)=-1^(i+k)/h0(k);
    end
end
Rx0=B0*Rjn*B0';
rx0d0=B0*Rjn*h0;
Wx0=inv(Rx0)*rx0d0;

%%%%%%%%%%%%%%MWF_SIM%%%%%%%%%%%%%%%
%前向递推
r=M-1;       % r MWF的阶数，三种算法结果完全一样，说明MWF与MVDR的等效性
h=[];
% Pd=[];
delta=[];
L=[];
% D=[];
B=eye(M-1);
rxd=rx0d0;
Rxi=Rx0;
for i=1:r
    if i<r
        deltai=(rxd'*rxd)^0.5;
        delta=[delta,deltai];
        hi=rxd/deltai;
        h=[h;hi];
        if i==1
            L=[L;hi'];
        else
            L=[L;hi'*B];
        end
        Bi=zeros(M-i-1,M-i);
        for kk=1:M-i-1
            for k=kk:kk+1
                Bi(kk,k)=-1^(kk+k)/hi(k);
            end
        end
        B=Bi*B;
        rxd=Bi*Rxi*hi;
%         Pdi=hi'*Rxi*hi;
%         Pd=[Pd,Pdi];
        Rxi=Bi*Rxi*Bi';
    else
        L=[L;B];
    end
end

Rd=L*Rx0*L';
rdd0=L*rx0d0;
Wd=inv(Rd)*rdd0;

[V,D]=eig(Rjn);
%对特征值从小到大进行排列并对对应的特征向量进行调整
for rr=1:M
for i=1:M-1
    if D(i,i)>D(i+1,i+1)
        Dtemp=D(i+1,i+1);
        D(i+1,i+1)=D(i,i);
        D(i,i)=Dtemp;
        Vtemp=V(:,i+1);
        V(:,i+1)=V(:,i);
        V(:,i)=Vtemp;
    end
end
end
Us=V(:,(M-numJam+1:M));
Ds=D((M-numJam+1:M),(M-numJam+1:M));
T=[Us,a0];

Un=V(:,(1:M-numJam));
Dn=D((1:M-numJam),(1:M-numJam));

dd=diag(D);
[dd,IX]=sort(dd);
dd=flipud(dd);
IX=flipud(IX);
U=[];
for i=1:r
    U=[U,V(:,IX(i))];
end


MVDR=1;
if MVDR==1
    W1=inv(Rjn)*a0/(a0'*inv(Rjn)*a0);
end


% ang=-90:0.5:90;
% a=zeros(M,1);
% for k=1:length(ang)
%     for i=1:M
%         a(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang(k)*pi/180));
%     end
%     beam1(k)=W1'*a;            
% end
% beam1=abs(beam1);
% beam1=beam1/max(beam1);
% beam1=20*log10(beam1);
% 
% figure(1)
% plot(ang,beam1);
% axis([-90  90  -50  0]);
% hold on
% plot([15 15], [-50 0], 'm','Linewidth',2);
% xlabel('Angle of Arrive(°)');
% hold off
% grid on


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%上行波束赋形直接用于下行

Mod = 8;                      % Modulation order
% hModem = modem.pskmod('M',8,'PhaseOffset', -pi/8,'SymbolOrder','Gray');   % 8-PSK modulator object
hModem = modem.pskmod(Mod);
hDeModem=modem.pskdemod(hModem);

Rsym = 9600;                % Input symbol rate
Rbit = Rsym * log2(M);      % Input bit rate
Nos = 4;                    % Oversampling factor
ts = (1/Rbit) / Nos;        % Input sample period

% fc=2e9;
fup = fc;                % Uplink Carrier frequency
freqGap=0;            %上下行频率间隔
fdown=2e9+freqGap;
c = 3e8;                    % Speed of light in free space
fs=5*fc;                    %采样频率
v = 0 * 1e3/3600;           % Mobile speed (m/s)
fd = v*fc/c;   


ArrayNum=4;                   %阵元数目
PathNum=1;                    %路径数目
lamdaD=c/fdown;

doa_theta=15*pi/180;       %入射方位角     
ang=sin(doa_theta);

nn=1:ArrayNum;
arrayvec=exp(j*2*pi/lamdaD*d*(nn-1)*sin(ang0*pi/180));
ArrayVec=arrayvec.';


Nburst = 4000;        %突发数
Nbits =20;            %1个突发150bit
Ts=ts*Nbits;

h1=rayleighchan(Ts,fd);
h1.StorePathGains = 1;          % Enables storage of path gains
h1.ResetBeforeFiltering = 0;    % The channel is not reset between frames
h1.NormalizePathGains = 1;    

SNR=0:1:10;
for n=1:length(SNR)
    noe=0;
    for iTxSig=1:Nburst
        h1.ResetBeforeFiltering = 0;   
        chan=h1.PathGains;
        fade=chan*chan';
        
        tdata=randint(Nbits,1,8);
        Tx=modulate(hModem,tdata);  %调制信号   
        x=Tx*W1';
        TxSig=x*ArrayVec;
        
%         fadeSig=filter(h1,TxSig);
%         SigPower=fadeSig'*fadeSig/Nbits;
%         randnoi=(sqrt(2)/2)*(randn(Nbits,1)+j*randn(Nbits,1));
%         noicor=SigPower/(10^(0.1*SNR(n)));
%         addNoi=noicor*randnoi;
%         
%         RxSig=(1/fade)*chan'*fadeSig;
%         rxSig=RxSig+addNoi;
        
        fadeSig=filter(h1,TxSig);
        
        %%%%%%%%%%%%%%%%%%%%%经过高斯信道
        SigPower=fadeSig'*fadeSig/Nbits;
        randnoi=(sqrt(2)/2)*(randn(Nbits,1)+j*randn(Nbits,1));
        noicor=SigPower/(10^(0.1*SNR(n)));
        addNoi=noicor*randnoi;
        RxSig=fadeSig+addNoi;
        
        %%%%%%%%%%%%%%%%%%%%%接收
        rxSig=(1/fade)*chan'*RxSig;
        rdata = demodulate(hDeModem, rxSig);
%         rdata = demodulate(hDeModem, TxSig);
        %%%%%%%%%%%%信道衰落参数变化
        h1.ResetBeforeFiltering = 1;
        y2=filter(h1,tdata);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%误码率计算
        for m=1:Nbits
            Wdata=0;
            if rdata(m)~=tdata(m)
                Wdata=1;
            end
            noe=noe+Wdata;
        end 
    end
     ser(n)=noe/(Nburst*Nbits);
end

figure(2)
semilogy(SNR,ser,'r-*');
hold on
% legend('FDD:frequency gap-45MHz');
% xlabel('SNR(dB)'); ylabel('SER');
% title('smart antennas in FDD over Rayleigh Fading Channel');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%FDD system,fu=2000MHz,fdown=2095MHz
M=4;             %array_number;
fc=2e9;          %work_frequency;
c=3e8;           %rate of light
lamda=c/fc;      %work_wavelength;
d=lamda/2;       %array_distance;

SNR=0;
JNR1=30;
JNR2=30;
JNR3=30;

ang0=15;
ang1=25;
ang2=0;
ang3=-15;

K=100;           %snap_number;
fs=1.8e10;       %sample_frequency;

a0=zeros(M,1);        % a0 期望信号方向矢量
a1=zeros(M,1);        % a1 干扰信号1的方向矢量
a2=zeros(M,1);        % a2 干扰信号2的方向矢量
a3=zeros(M,1);        % a3 干扰信号3的方向矢量
j=sqrt(-1);
for i=1:M
    a0(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang0*pi/180));
    a1(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang1*pi/180));
    a2(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang2*pi/180));
    a3(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang3*pi/180));
end
% A=[a0';a1';a2';a3']';
% A1=[a1';a2';a3']';
A=[a0,a1,a2,a3];
A1=[a1,a2,a3];

t=(1:1:K)/fs; 

s0=cos(2*pi*fc*t);        

AM_JAM1=1;
PM_JAM1=0;
AM_JAM2=1;
PM_JAM2=0;
AM_JAM3=1;
PM_JAM3=0;


if AM_JAM1==1
    s01=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM1==1
    s01=sin(2*pi*fc*t+randn(1,K));
end

if AM_JAM2==1
    s02=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM2==1
    s02=sin(2*pi*fc*t+randn(1,K));
end

if AM_JAM3==1
    s03=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM3==1
    s03=sin(2*pi*fc*t+randn(1,K));
end


Ps0=s0*s0'/K;                       % 计算期望信号的功率
Ps1=s01*s01'/K;                     % 计算干扰信号1的功率
Ps2=s02*s02'/K;                     % 计算干扰信号2的功率
Ps3=s03*s03'/K;                     % 计算干扰信号2的功率

Ps=10.^(SNR/10);
Pj1=10.^(JNR1/10);
Pj2=10.^(JNR2/10);
Pj3=10.^(JNR3/10);

s=s0/Ps0^0.5*Ps^0.5;             % 具有输入信噪比的期望信号
s1=s01/Ps1^0.5*Pj1^0.5;          % 具有干噪比1的干扰信号1
s2=s02/Ps2^0.5*Pj2^0.5;          % 具有干噪比2的干扰信号2
s3=s03/Ps3^0.5*Pj3^0.5;          % 具有干噪比2的干扰信号2

S=[s;s1;s2;s3];
S1=[s1;s2;s3];
numJam=size(S1,1);
Noise=randn(M,K);
X=A*S+Noise;
X1=A1*S1+Noise;

R=X*X'/K;
Rjn=X1*X1'/K;


h0=a0/(a0'*a0)^0.5;
B0=zeros(M-1,M);
for i=1:M-1
    for k=i:i+1
        B0(i,k)=-1^(i+k)/h0(k);
    end
end
Rx0=B0*Rjn*B0';
rx0d0=B0*Rjn*h0;
Wx0=inv(Rx0)*rx0d0;

%%%%%%%%%%%%%%MWF_SIM%%%%%%%%%%%%%%%
%前向递推
r=M-1;       % r MWF的阶数，三种算法结果完全一样，说明MWF与MVDR的等效性
h=[];
% Pd=[];
delta=[];
L=[];
% D=[];
B=eye(M-1);
rxd=rx0d0;
Rxi=Rx0;
for i=1:r
    if i<r
        deltai=(rxd'*rxd)^0.5;
        delta=[delta,deltai];
        hi=rxd/deltai;
        h=[h;hi];
        if i==1
            L=[L;hi'];
        else
            L=[L;hi'*B];
        end
        Bi=zeros(M-i-1,M-i);