all.m

OFDM频偏估计优化算法代码

%
%
%研学论坛上的程序，用三种不同的方法（MS/MMSE/ML）估计cp同步的位置

profile on 

% Initialize the parameters
NumLoop = 10000;
N = 256; %the number of the subcrarries
L = 32; % the length of CP
ErrSym=zeros(3,length(0:2:20));
count_err_matrix = zeros(3,NumLoop,length(0:2:20));%生成三维矩阵，用来装载fft估计窗口起始值与真实值得差值
                                                   %行装载三种不同的方法，列装载每个ofdm符号估计的误差，径装载不同的信噪比下的误差

for SNR=0:2:20%loop and calculate differint SNR
    temp=[0 0 0];
    SNRtmp=SNR/(1+SNR);
    
     for loop=1:NumLoop%main loop
        
% Generate the random binary stream for transmit test
BitsTx = floor(rand(1,N*4)*2);%1024 bits

% Modulate (generates QAM symbols)
SymQAMtmp = reshape(BitsTx,2,N*2).';%finally the matrix is N*2 by 2
%RESHAPE(X,M,N) returns the M-by-N matrix whose elements
%    are taken columnwise from X.  An error results if X does
%    not have M*N elements.


SymQAMtmptmp = bi2de(SymQAMtmp,2,'left-msb');%the base is 2 , not 10
%D = BI2DE(B) converts a binary vector B to a decimal value D. When B is a
%  matrix, the conversion is performed row-wise and the output D is a column
%    vector of decimal values. The default orientation of the binary input
%    is Right-MSB; the first element in B represents the least significant bit.


%00->-1-i,01->-1+i,10->1-i,11->1+i
QAMTable = [-1-i -1+i 1-i 1+i];
SymQAM = QAMTable(SymQAMtmptmp+1);%genetate QAM symbels

% Add cyclic prefix
NumAddPrefix = 2*N + L;
SymCP = zeros(1,NumAddPrefix);
RowPrefix = (2*N - L + 1):2*N;
SymCP = [SymQAM(:,1:N/2),SymQAM(:,RowPrefix),SymQAM(:,N+1:2*N),SymQAM(:,N/2+1:N)];%commen trasmite OFDM symbels with one cp in middel position
%set the lengh of the data is 2N+L

% Go through the channel
% input: SymCP(N + L,NumLoop); output: SymCh(1,(N + L)*NumLoop)
SymCh = zeros(1,2*N + L);
SymChtmp = SymCP;
Ch = [1 1/2 1/4];
SymChtmptmp = filter(Ch,1,SymChtmp);%2*N+L lenght
% Add the AWGN
SymCh = awgn(SymChtmptmp,SNR,'measured');
%Y = AWGN(X,SNR,SIGPOWER) when SIGPOWER is numeric, it represents 
%    the signal power in dBW. When SIGPOWER is 'measured', AWGN measures
%    the signal power before adding noise.
RxSignal=SymCh;%2*N+L lenght
S=zeros(1,(length(RxSignal)-N-L+1));
P=zeros(1,(length(RxSignal)));
ML=S;
MMSE=ML;
for k=1:(length(RxSignal)-N-L+1)%one ofdm symble(uninclude cp) lenght  N
    rtmp=RxSignal(1,k:k+L-1);%32 chips lenght slide wendow
    rtmptmp=RxSignal(1,k+N:k+L+N-1);%N chips delay wendow
  % rtmptmp=conj(rtmptmp);
    S(k)=rtmp*rtmptmp';%mutiplay and sum
    P(k)=rtmp*rtmp';
    S(k)=abs(S(k));
end
for k=1:(length(RxSignal)-N-L+1)%N times loop
   MMSE(k)=S(k)-1/2*(P(k+N)+P(k));
   ML(k)=S(k)-1/2*SNRtmp*(P(k+N)+P(k));
end
[numtmp11 numtmp12]=max(S);
count_err_matrix(1,loop,(SNR/2 + 1)) = abs(numtmp12 - (N/2 + 1));%calculate error 
[numtmp21 numtmp22]=max(MMSE);
count_err_matrix(2,loop,(SNR/2 + 1)) = abs(numtmp22 - (N/2 + 1));
[numtmp31 numtmp32]=max(ML);
count_err_matrix(3,loop,(SNR/2 + 1)) = abs(numtmp32 - (N/2 + 1));


if numtmp12 ~= N/2+1
   temp(1)=temp(1)+1;
end
if numtmp22 ~= N/2+1
   temp(2)=temp(2)+1;
    end
if numtmp32 ~= N/2+1
  temp(3)=temp(3)+1;
end
end
ErrSym(:,SNR/2+1)=[temp(1)/NumLoop;temp(2)/NumLoop;temp(3)/NumLoop];
end
figure(1),plot(0:2:20,ErrSym(1,:),'x:',0:2:20,ErrSym(2,:),'r*:',0:2:20,ErrSym(3,:),'bo:');
xlabel('SNR(db)');ylabel('符号同步错误率');
title('OFDM符号同步');
legend('MC','MMSE','ML');


%function [ output_args ] = Untitled1( input_args )
%UNTITLED1 Summary of this function goes here
%  本程序对myofdm程序中产生的fft窗口估计值与实际值的差值矩阵进行统计处理。
%计算出差值矩阵的在不同方法（）不同信噪比（）下的偏差点数目的分布
%计算出偏差点的统计方差矩阵和均值矩阵，并绘图显示。


%load dbb_myofdm.mat;

%计算错误计数矩阵的错误估计值的相应的错误数目
%
%

MS_count_err_matrix = reshape(count_err_matrix(1,:,:),NumLoop,length(0:2:20));%三维矩阵转二维矩阵
MMSE_count_err_matrix = reshape(count_err_matrix(2,:,:),NumLoop,length(0:2:20));
ML_count_err_matrix = reshape(count_err_matrix(3,:,:),NumLoop,length(0:2:20));

MS_count_err_matrix = MS_count_err_matrix';%将各算法的计数误差矩阵转换为与三维矩阵统一的格式
MMSE_count_err_matrix = MMSE_count_err_matrix';
ML_count_err_matrix = ML_count_err_matrix';

MS_sort = sort(MS_count_err_matrix,2);%二维矩阵按升序排序
MMSE_sort = sort(MMSE_count_err_matrix,2);
ML_sort = sort(ML_count_err_matrix,2);

MS_max_num = max(MS_sort(:,NumLoop));%矩阵中的最大值
MMSE_max_num = max(MMSE_sort(:,NumLoop));
ML_max_num = max(ML_sort(:,NumLoop));

MS_statistic_matrix = zeros( length(0:2:20),(MS_max_num+1) );%生成MS统计矩阵规模
for m = 1:length(0:2:20);%不同信噪比下的数值
    for n = 0:1:MS_max_num;%不同差值出现的数目求和统计
        index = find(MS_sort(m,:) == n);
        total = length(index);
        MS_statistic_matrix(m,(n+1)) = total;
    end
end

MMSE_statistic_matrix = zeros( length(0:2:20),(MMSE_max_num+1) );%生成MMSE统计矩阵规模
for m = 1:length(0:2:20);%
    for n = 0:1:MMSE_max_num;%
        index = find(MMSE_sort(m,:) == n);
        total = length(index);
        MMSE_statistic_matrix(m,(n+1)) = total;
    end
end

ML_statistic_matrix = zeros( length(0:2:20),(ML_max_num+1) );%生成ML统计矩阵规模
for m = 1:length(0:2:20);%
    for n = 0:1:ML_max_num;%
        index = find(ML_sort(m,:) == n);
        total = length(index);
        ML_statistic_matrix(m,(n+1)) = total;
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%计算不同信噪比下的错误估计的方差
%
%
MS_snr_standard = std((MS_count_err_matrix+(N/2+1)),0,2);%
MMSE_snr_standard = std((MMSE_count_err_matrix+(N/2+1)),0,2);
ML_snr_standard = std((ML_count_err_matrix+(N/2+1)),0,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算不同信噪比下的错误估计的均值
%
%
MS_snr_mean = mean(MS_count_err_matrix,2);%
MMSE_snr_mean = mean(MMSE_count_err_matrix,2);
ML_snr_mean = mean(ML_count_err_matrix,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%画图输出
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)%MS方法的偏差点数目的分布
stem3(MS_statistic_matrix);
axis([0 130 0 12 -inf inf]);
title('MS');
xlabel('error numerical value of estimation');
ylabel('differint SNR');
zlabel('total numbers');

figure(3)%MMSE方法的偏差点数目的分布
stem3(MMSE_statistic_matrix);
axis([0 130 0 12 -inf inf]);
title('MMSE');
xlabel('error numerical value of estimation');
ylabel('differint SNR');
zlabel('total numbers');

figure(4)%ML方法的偏差点数目的分布
stem3(ML_statistic_matrix);
axis([0 130 0 12 -inf inf]);
title('ML');
xlabel('error numerical value of estimation');
ylabel('differint SNR');
zlabel('total numbers');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(5)%不同信噪比下的不同方法的方差
plot(MS_snr_standard,'-ro');
hold on
plot(MMSE_snr_standard,'--b+');
plot(ML_snr_standard,':k*');
hold off
xlabel('SNR');
ylabel('error standard');
title('MS/MMSE/ML error standard');
legend('MS','MMSE','ML');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(6)%不同信噪比下的不同方法的均值
plot(MS_snr_mean,'-ro');
hold on
plot(MMSE_snr_mean,'--b+');
plot(ML_snr_mean,':k*');
hold off
xlabel('SNR');
ylabel('error mean');
title('MS/MMSE/ML error mean');
legend('MS','MMSE','ML');






profile off
profile viewer
profsave(profile('info'),'profile_results_myofdm')