estcorr.m

基于卷积信号的MIMO系统盲信号估计

function [pow_est, pow_true]=estCorr(NF, H, h, counter,SNR,length_diff)

%m: # of sensors
%n: # of sources
[L,m,n]=size(h);

if mod(NF,2) KF=(NF+1)/2; else KF=NF/2; end

N=4048*4;                       %%% Length of Input(Output) signal 
seg_length = 4048*2;            %%% Segment length used in estimating the cross cumulants
seg_num=N/seg_length;            %%% Number of segments
C_LENGTH=L+length_diff;                    %%% how much cumulants need estimate, maximum argument of cumulants

for ii=1:NF
    pow_true(ii,:,:)=H(:,:,ii)*H(:,:,ii)';
end


%%%Estimation Part


rand('seed',counter)
s = -log(rand(n,N));      %%% Single side Exponential distributed real signal
s = s-(mean(s'))'*ones(1,N);
s=diag(1./std(s,0,2))*s;       %% Input Signal
x = zeros(size(s));          %%% Observed mixture , received signal
for ii=1:m
    for jj=1:n
        x(ii,:) = x(ii,:) + filter(h(:,ii,jj),1,s(jj,:));
    end               
end
x = x-(mean(x'))'*ones(1,N);

for ii=1:m

    std_x(ii)=std(x(ii,:));
    
end

    
for ii=1:m
        
    rand('state',sum(100*clock)*rand);
    noise=randn(1,N);
    noise=noise-mean(noise);
    power_coefficient = (10^(SNR/20))*std(noise)/std(x(ii,:));
    x(ii,:) = x(ii,:)+noise/power_coefficient;
    
    x(ii,:)=std_x(ii)*x(ii,:)/std(s(ii,:));
end

corr_est = zeros(2*C_LENGTH+1,m,m);

for (ii=1:m)
        for (jj=1:m)
                for seg=0:seg_num-1
                        
                    tmp=xcorr(x(ii,seg_length*seg+1:seg_length*(seg+1)),...
                            x(jj,seg_length*seg+1:seg_length*(seg+1)),C_LENGTH,'unbiased')/seg_num;
                    
                    tmp=reshape(tmp,[2*C_LENGTH+1,1,1]);
                        
                    corr_est(:,ii,jj)=corr_est(:,ii,jj)+tmp;
                 end
          end
 end


pow_est=zeros(NF,m,m);
for ii=1:m
    for jj=1:m
            cx_dummy=zeros(NF,1);
            cx_dummy((KF+1-C_LENGTH) : (KF+1+C_LENGTH))=corr_est(:,ii,jj);
            cx_dummy=fftshift(cx_dummy);
            B_x_dummy=(fft(cx_dummy,NF));
            pow_est(:,ii,jj)=B_x_dummy;
    end
end

       if 0
            figure; clf

             for ii=1:m
                   for jj=1:m
                           ratio=max(abs(squeeze(pow_true(:,ii,jj))))/max(abs(squeeze(pow_est(:,ii,jj))));
                           subplot(m,m,jj+(ii-1)*m), plot(abs(squeeze(pow_est(:,ii,jj)))*ratio); hold on;
                           plot(abs(squeeze(pow_true(:,ii,jj))),'r'); hold on;
                    end
                end
            
            
        end