responeofchannel.m

此程序是关于OFDM（正交频分复用）同步中信道响应的程序

%信道估计。可以发现，频率偏移被纠正后得到的信道的估计，和只有噪声没有频偏时的信道估计是一样的。没有频偏时程序带着
%频偏估计的过程对结果没有影响。
%    利用方法二取信道的平均时，注意将频偏设为零进行。这时相邻符号的初始相位同为0，
%才能取平均。  这是为了编程方便，事实上在进行频偏补偿时如果对两个符号一起补偿，则其初始相位相同，就可以进行平均。
%

clear all;
EstimateMethod=1;    %只利用第二个符号估计信道
EstimateMethod=2;    %利用两个符号分别估计信道，并进行平均。
hn = zeros(1,31);
hn(1) = (1+j);
hn(8) = (0.5-0.3*j);
hn(10) = (0.3+0.2*j);
hn(22) = (0.2-0.1*j);
hn(31) = 0.1-0.1j;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters can be adjusted
NumofSubcarrier = 256;    
CyclicExtentionLength = 32;
NumOfSymbol = 2000;                     
MappingScheme = 'BPSK';    %这里映射方式没有实际的用处，只是为了调用通用的信道函数。频域的伪随机序列正好相当于BPSK。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PNCode = GenPNCode(NumofSubcarrier);
PNCodeTimeDomain = ifft(PNCode);
PNCodeTimeDomain = PNCodeTimeDomain./sqrt(sum(abs(PNCodeTimeDomain).^2)/NumofSubcarrier);   %能量归一化
TransData = [];
for n=1:1:NumOfSymbol
  TransData = [TransData,PNCodeTimeDomain.'];   %直接构造时域符号形式，每列是一个符号
end
ParallelSample = TransData;  
ParallelSample = AddCyclicExtention(ParallelSample,CyclicExtentionLength);  %如果使用循环前缀，则加此句，否则去掉
SerialSample = ParallelToSerial(ParallelSample);                                        %

LogEbNo = 30;
FrequencyOffset = 0;

ChannelEstimate = zeros(NumOfSymbol-1,NumofSubcarrier);

RecieveData = MultipathChannel(SerialSample,LogEbNo,hn,NumOfSymbol,MappingScheme);               %pass througn a MultipathChannel

DataWithOffset = AddContinusOffset(RecieveData,NumofSubcarrier,FrequencyOffset);

for q=0:1:NumOfSymbol-2
   A = DataWithOffset(q*(NumofSubcarrier+CyclicExtentionLength)+CyclicExtentionLength +2:(q+1)*(NumofSubcarrier+CyclicExtentionLength)).*conj(DataWithOffset((q+1)*(NumofSubcarrier+CyclicExtentionLength)+CyclicExtentionLength +2:(q+2)*(NumofSubcarrier+CyclicExtentionLength)));
   A_angle = atan(-1*sum(imag(A))/sum(real(A)));
   if sum(real(A)) < 0
     A_angle = A_angle + pi;
     if A_angle > pi
        A_angle = A_angle - 2*pi;
     end
   end
   OffsetEstimate = A_angle*NumofSubcarrier/(2*pi*(NumofSubcarrier+CyclicExtentionLength));
   DataWithouOffset = AddContinusOffset(DataWithOffset(q*(NumofSubcarrier+CyclicExtentionLength)+1:(q+1)*(NumofSubcarrier+CyclicExtentionLength)),NumofSubcarrier,-1*OffsetEstimate);
   ChannelEstimate(q+1,:) = fft(DataWithouOffset(CyclicExtentionLength+1:end)).*PNCode;
end

VarChannelEstimate = zeros(1,NumofSubcarrier);
Hn = fft([hn,zeros(1,NumofSubcarrier-length(hn))]);
Hn = Hn./sqrt(sum(abs(Hn).^2)./NumofSubcarrier);         %能量归一化
Temp1=0;
Temp2=0;
for m=1:1:NumOfSymbol-1
   ChannelEstimate(m,:) = ChannelEstimate(m,:)./sqrt(sum(abs(ChannelEstimate(m,:)).^2)/NumofSubcarrier);  %能量归一化
   if EstimateMethod == 1
      VarChannelEstimate = VarChannelEstimate + (abs(ChannelEstimate(m,:))-abs(Hn)).^2;
   else
      Temp2 = ChannelEstimate(m,:);
      if m > 1   %第一个值由于无法和前面的平均而出现异常，不用。
        VarChannelEstimate = VarChannelEstimate + (abs((Temp2+Temp1)./2)-abs(Hn)).^2;
      end
      Temp1=Temp2;
   end
end

if EstimateMethod == 1
   VarChannelEstimate = VarChannelEstimate/(NumOfSymbol-1);
else
   VarChannelEstimate = VarChannelEstimate/(NumOfSymbol-2);   %比上面的算法少了第一个无法计算的异常值。
end

plot(VarChannelEstimate,'r');


%end of the program