blocks = 200;      %OFDM塊數
%與chanestnew相比，是將線性相關改為圓周相關
%梳狀導頻與本文方法（迭代前和迭代后）作比較
% |--cp--|----PN----|------OFDM------|
tic;
fc = 2000;       %MHz,系統載波頻率
N = 512;   %子載波數
L = 32;    %循環前綴長度
Nf3 = 6;    %導頻率為
M = 26;     %多徑數%%%%%%%%%%%%%%

%generate transmitted signal   s   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DN = blocks*N;
for j = 1:DN
    temp = rand;
    if (temp<0.25)
        s(j) = 1;
    elseif (temp<0.5)
        s(j) = sqrt(-1);
    elseif (temp<0.75)
        s(j) = -1;
    else 
        s(j) = (-1)*sqrt(-1);
    end
end

%generate pn code    pns   1*pL   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pn = genpncode(6);
pL = length(pn)+L;
pns = [pn pn(1:L)];
PNS = fft(pns,N)/sqrt(N);

%transmit & receive    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n = 1:blocks                      
    for k = 1:N
        TR(n,k) = s(N*(n-1)+k);       %convert series to parallel
    end 
end
for n = 1:blocks
    tr1(n,:) = [pns ifft(TR(n,pL+1:N))*sqrt(N-pL)];
end
tr1 = [tr1(:,(N-L+1):N) tr1];       %blocks*(N+L)
for j = 1:blocks
    ss1((j-1)*(N-pL)+1:j*(N-pL)) = s((j-1)*N+pL+1:j*N);                 %本文所提方法的傳輸信號
end

pfc = sqrt(-1)*ones(blocks,1);    %頻域梳狀導頻方法3發射的數據 
TR3 = TR;
for j = 1:Nf3:N-2%%%%%%%%%%%%%
    TR3(:,j) = pfc;
end
TR3(:,end) = pfc;%%%%%%%%%%%%%%%%%%%%%%%%%%
for n = 1:blocks
    tr3(n,:) = sqrt(N)*ifft(TR3(n,:))
end
tr3 = [tr3(:,(N-L+1):N) tr3];       %blocks*(N+L)
j = 1;
Npc = round(N/Nf3)+1;      %梳狀導頻數 Npc = (N-1)/Nf3+1;%%%%%%%%%%%%%%
%for k = 1:N
    %if (mod(k-1,Nf3)~=0)
       %ss3(j:(Npc-1)*(Nf3-1):j+(blocks-1)*(Npc-1)*(Nf3-1)) = TR3(:,k);  %得到除導頻信號外的數據信號
       %j = j+1;
    %end
%end
for k = 1:N-1%
    if (mod(k-1,Nf3)~=0)%
        ss3(j:(Npc-1)*(Nf3-1)+1:j+(blocks-1)*((Npc-1)*(Nf3-1)+1)) = TR3(:,k);   %%得到除導頻信號外的數據信號
        j = j+1;%
    end%
end%
ss3(j:(Npc-1)*(Nf3-1)+1:j+(blocks-1)*((Npc-1)*(Nf3-1)+1)) = TR3(:,N-1); %

%delayTime = [0 0.2 0.6 1.6 2.4 5.0]*10^(-6);%%%%%%%%%%%%%%%%%%%
%delayn = delayTime/chipT;%%%%%%%%%%%%%%%%%%%%%
%averPower = %%%%%%%%%%%%%%%%%%

%channel %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
pnss = pns(L+1:pL);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delayn = 0:25;%[0 1 2 3];%[0 1 3 8 12 25];
%averPower = [0.189 0.379 0.239 0.095 0.060 0.038];[0.575 0.362 0.057 0.006];
snrdB = 0:5:25;
for m = 1:length(snrdB)
    %%%%%%%%%%%%%%%%%%%%%%%channel model
    h = zeros(N,blocks);%
    for k=1:M%
        h(delayn(k)+1,:)=randn(1,blocks)+sqrt(-1)*randn(1,blocks);%
    end%
    
    for k = 1:M%
        h(delayn(k)+1,:) = h(delayn(k)+1,:)/norm(h(delayn(k)+1,:))*sqrt(blocks);%             %歸一化信道 N*blocks
    end%
    %for k = 1:M
        %h(delayn(k)+1,:) = sqrt(averPower(k))*h(delayn(k)+1,:);
    %end
    for n = 1:blocks
        h(:,n) = h(:,n)/norm(h(:,n));%歸一化信道
    end
    HL = fft(h);    %  N*blocks
    H(:,:,m) = HL';
    Hs = reshape(H(:,:,m),1,N*blocks);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
    h1 = zeros(blocks,N);
    h2 = zeros(blocks,N);%參數檢測未改進之前
    for n = 1:blocks
        r1(n,:) = zeros(1,N);
        r3(n,:) = zeros(1,N);
        for j = 1:M
            r1(n,:) = r1(n,:)+h(j,n).*tr1(n,(L+1-(j-1):N+L-(j-1)));%
            r3(n,:) = r3(n,:)+h(j,n).*tr3(n,(L+1-(j-1):N+L-(j-1)));%
        end
        r1(n,:) = AWGN(r1(n,:),snrdB(m));
        r3(n,:) = AWGN(r3(n,:),snrdB(m));
        R1(n,:) = fft(r1(n,:))/sqrt(N);        %blocks*N
        %R12(n,:) = fft(r12(n,:))/sqrt(N);
        R3(n,:) = fft(r3(n,:))/sqrt(N);        %blocks*N
        
        Rprw = circlecorr(pnss,r1(n,L+1:pL));   %為本文方法1作準備
        Rs = Rprw;
        for j = 1:M;
            Rpp = Rs(pL-L:-1:1);
            [a(n,j,m),b(n,j,m)] = max(Rpp);
            a(n,j,m) = a(n,j,m)/(pL-L);
            pnsd = [pnss(pL-L-b(n,j,m)+2:pL-L) pnss(1:pL-L-b(n,j,m)+1)];
            Rppd = circlecorr(pnss,pnsd);
            Rs = Rs-a(n,j,m)*Rppd;
            h1(n,b(n,j,m)) = a(n,j,m);    %本文所提方法結果
        end
        H1(n,:) = fft(h1(n,:))/sqrt(N);
            
        Rs2 = Rprw;   %為未迭代方法作準備
        for j = 1:M;
            Rpp2 = Rs2(pL-L:-1:1);
            [c(n,j,m),d(n,j,m)] = max(Rpp2);
            c(n,j,m) = c(n,j,m)/(pL-L);
            Rs2(pL-L-d(n,j,m)+1) = 0;
            h2(n,d(n,j,m)) = c(n,j,m);    %本文所提方法結果
        end
        H2(n,:) = fft(h2(n,:))/sqrt(N);
    end
    
    k3 = 1;
    for j = 1:N-2;%j = 1:N
        if (mod(j-1,Nf3)==0)
           H3(:,k3) = R3(:,j)./pfc;         %頻域梳狀導頻估計3結果
           k3 = k3+1;
        end
    end
    H3(:,k3) = R3(:,N)./pfc; %%%%%%%%%%%%%%%%%  
    
    [X3,Y3] = meshgrid([1:Nf3:N-2 N],1:blocks);%[X3,Y3] = meshgrid(1:Nf3:N,1:blocks);%%%%%%%%%%%
    [X,Y] = meshgrid(1:N,1:blocks);
    HH1(:,:,m) = H1*sqrt(N);    %本文方法1 插值完畢blocks×N
    HH2(:,:,m) = H2*sqrt(N);
    HH3(:,:,m) = interp2(X3,Y3,H3,X,Y,'linear');    %頻域梳狀導頻方法4插值完畢

    %均衡＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
    TRe1 = R1./HH1(:,:,m);
    for n = 1:blocks
        tre1(n,:) = ifft(TRe1(n,:)-PNS)*sqrt(N);
        TRe1(n,:) = [zeros(1,pL) fft(tre1(n,pL+1:N))/sqrt(N-pL)];
        TRes1(N*(n-1)+1:N*n) = TRe1(n,:);        %并串轉換
    end
    for j = 1:blocks
        TRess1((j-1)*(N-pL)+1:j*(N-pL)) = TRes1((j-1)*N+pL+1:j*N);                 %本文所提方法均衡后的頻域傳輸信號
    end
    
    TRe2 = R1./HH2(:,:,m);
    for n = 1:blocks
        tre2(n,:) = ifft(TRe2(n,:)-PNS)*sqrt(N);
        TRe2(n,:) = [zeros(1,pL) fft(tre2(n,pL+1:N))/sqrt(N-pL)];
        TRes2(N*(n-1)+1:N*n) = TRe2(n,:);        %并串轉換
    end
    for j = 1:blocks
        TRess2((j-1)*(N-pL)+1:j*(N-pL)) = TRes2((j-1)*N+pL+1:j*N);                 %迭代前本文所提方法均衡后的頻域傳輸信號
    end
    
    j = 1;
    TRe3 = R3./HH3(:,:,m);
    %for k = 1:N
        %if (mod(k-1,Nf3)~=0)
            %TRess3(j:(Npc-1)*(Nf3-1):j+(blocks-1)*(Npc-1)*(Nf3-1)) = TRe3(:,k);   %
            %j = j+1;
        %end
    %end
    for k = 1:N-1
        if (mod(k-1,Nf3)~=0)
            TRess3(j:(Npc-1)*(Nf3-1)+1:j+(blocks-1)*((Npc-1)*(Nf3-1)+1)) = TRe3(:,k);   %
            j = j+1;
        end
    end
    TRess3(j:(Npc-1)*(Nf3-1)+1:j+(blocks-1)*((Npc-1)*(Nf3-1)+1)) = TRe3(:,N-1);
    
    %解調并計算誤比特率誤符號率
    [ps1(m),pb1(m)] = demanderr(ss1,TRess1);
    [ps3(m),pb3(m)] = demanderr(ss3,TRess3);
    [ps2(m),pb2(m)] = demanderr(ss1,TRess2);%
    
    HH1s = reshape(HH1(:,:,m),1,N*blocks);
    HH2s = reshape(HH2(:,:,m),1,N*blocks);
    HH3s = reshape(HH3(:,:,m),1,N*blocks);
    cc1(:,:,m) = corrcoef(abs(HH1s),abs(Hs));
    cc2(:,:,m) = corrcoef(abs(HH2s),abs(Hs));
    cc3(:,:,m) = corrcoef(abs(HH3s),abs(Hs));
end
simulation_T=toc;

figure;
semilogy(snrdB,pb2,'b-*',snrdB,pb1,'b-d',snrdB,pb3,'b-+');
legend('時域估計方法','時域估計方法（改進信道參數檢測）','頻域梳狀導頻方法');
xlabel('信噪比/dB');
ylabel('誤碼率');
grid on;
axis([0,25,10^(-4),10^0]);
figure;%
plot(1:N,abs(H(1,:,4)),'b-',1:N,abs(HH2(1,:,4)),'b--',1:N,abs(HH1(1,:,4)),'b:',1:N,abs(HH3(1,:,4)),'b-.');%
legend('真實信道','時域估計方法','時域估計方法（改進信道參數檢測）','頻域梳狀導頻方法');%
xlabel('子載波');%
ylabel('信道幅頻響應');%
axis([1 N 0 2]);