svm_qpsk.m

基于支持向量机的盲均衡算法,使用的输入信号是QPSK.

%常模信号支持向量机盲均衡算法
%信号为QPSK，参考文献I.Santanmaria "Blind Equalization of Constant Modulus 
%Signals Using Support Vector Machines" IEEE trans. on Signal Processing
%2004,52(6):1773-1781
clear all;
clc;
close all;
%%%%%%%%%%%%信道的选定
j=sqrt(-1);
% b=[0.7 -1];
% a=[1 -0.7];
% h=filter(b,a,[1,zeros(1,19)])*exp(j*pi/4);
% 
% h=[1 0.1294-0.483*j];
h=[0.0410+0.0109j 0.0495+0.0123j 0.0672+0.0170j 0.0919+0.0235j 0.7920+0.1281j 0.3960+0.0871j 0.2715+0.0498j 0.2291+0.0414j 0.1287+0.0154j 0.1032+0.0119j];
%------------参数设置---------%
Num_of_realizations=1; %学习次数
N=100;                  %信号长度；
SNR=30;                %输入信噪比
niter=30;              %IRWPQ最多迭代次数
K=31;                  %均衡器长度
%---支持向量回归（SVR的参数------%
lamda=0.9;            %式22的参数λ
loss='eInsensitive';  %Vapnik ξ不敏感损失函数
%loss='quadratic';
ker='linear';         %线性内核
ISI_tot=[];
AME_tot=[];           

for t=1:Num_of_realizations
    t
    isi=[];
    ame=[];
    s=qpsk(N+K);
    x1=filter(h,1,s(:));
    noise_var=10^(-SNR/10)*sum(abs(h).^2);
    noise=sqrt(noise_var/2)*randn(N+K,1)+j*sqrt(noise_var/2)*randn(N+K,1);
    x=x1+noise;
    figure(1)
    plot(real(x),imag(x),'*');%均衡器输入信号星座图
    title('均衡前的星座图');
    %-----------------------------SVR最优参数-------%
    C=mean(abs(x).^2)+3*std(abs(x).^2);%参数C
    e=3*sqrt(noise_var)*sqrt(log(N)/N);%参数ξ，τ=3
    %---------------------------%
    X=zeros(K,N);
    n=K+1;
    for f=1:N
        X(:,f)=x(n+f-1:-1:n+f-K);
    end
    %--------------------初始化---------------------%
    w=zeros(K,1);
    w((K+1)/2)=1;               %均衡器初始值
    beta=pinv(conj(X))*w;       %式22
    y=(conj(X)')*w;             %均衡器输出
    %---------------计算AME和ISI 初始值------------%
    aux1=conv(h,w);
    s_aux1=sum(abs(aux1).^2);
    m_aux1=max(abs(aux1).^2);
    isi_ini=10*log10((s_aux1-m_aux1)/m_aux1);
    isi=[isi isi_ini];
    ame=[ame sqrt(mean((abs(y)-1).^2))];
    %-------------------------------------
    true=1;
    iter=0;
    while(true==1)
        yr=real(y);yi=imag(y);
        Xr=real(X)';Xi=imag(X)';
        upper_part=Xr.*(yr*ones(1,K))+Xi.*(yi*ones(1,K));
        lower_part=-Xi.*(yr*ones(1,K))+Xr.*(yi*ones(1,K));
        G=[upper_part lower_part];      %%计算式11
        Y=ones(N,1);
        [nsv,alpha,bias]=svr_complex(G,Y,ker,C,loss,e);%计算式14
        beta_QP=alpha.*y;
        beta=lamda*beta+(1-lamda)*beta_QP;%式22
        w=conj(X)*beta;
        y=conj(X')*w;
        aux1=conv(h,w);
        s_aux1=sum(abs(aux1).^2);
        m_aux1=max(abs(aux1).^2);
        isi_fin=10*log10((s_aux1-m_aux1)/m_aux1);
        isi=[isi isi_fin];
        ame=[ame sqrt(mean((abs(y)-1).^2))];
        iter=iter+1;
        %----------------------退出循环标准------%
        if ((ame(end)-ame(end-1))>-1e-5)|(iter>50)
            true=0;
        end
    end
    figure(2)
     plot(real(y),imag(y),'*');%输出星座图
     title('均衡后星座图');

     figure(3)
     plot(isi);%剩余ISI收敛曲线
     ylabel('剩余ISI');
     xlabel('IRWQP迭代次数');
     figure(4)
    plot(ame);%平均模误差收敛曲线AME 见表达式24
    ylabel('AME');
     xlabel('IRWQP迭代次数');
end