OFDM 에서 CP 를 이용한 Time 과 Frequency Offset 검출 (MATLAB)

참고 논문 :

https://www.researchgate.net/publication/3316488_ML_Estimation_of_time_and_frequency_offset_in_OFDM_systems

 

CP 를 이용해 Time Offset 과 Freq Offset 을 검출하는 수식을

1. 직관적인 유도 (쉬움)

2. 가우시안 분포와 log-likelihood(log+조건부 확률)을 이용한 수식적 유도 (어려움)

두 과정으로 보겠습니다.

 

직관적인 유도

Time Offset

 

한 심볼의 길이가 N

CP 의 길이가 L 일 때

임의의 신호 부분에서 2N+L 만큼 추출 하면

반드시 신호 (N) 과 CP (L) 이 포함되어 있습니다.

 

CP 는 N 이후의 신호(CP' 라고 칭하겠음)에 대해 동일한 값을 가지고

신호는 복소 값들로 이루어져 있으므로

임의의 신호와 N 이후의 신호를 곱했을 때

CP 에 해당하는 부분이 최대 값을 가지게 됩니다.

 

또한 이 과정을 CP 의 길이 동안 하면

CP 의 시작점에서 곱한 결과를 더한 값에서

최대 값을 나타나므로

CP 의 시작 지점 또한 알 수 있습니다.

 

이 값은 Time Offset 을 의미합니다.

 

Frequency Offset

ofdm 에서 각 심볼은 주파수 영역에서 생성한 데이터를

IDFT(IFFT) 를 거쳐 시간 영역으로 전환해 데이터를 전송합니다.

IDFT 의 식은 아래와 같습니다.

자연상수 e 의 승수가 위상의 회전을 의미 하는 것을 생각하면

IDFT 를 통해 나오는 수열의 각 요소들은

원에서 일정한 간격(2 pi / N)을 두고 배치 된 것임을 알 수있습니다.

또한 해당 신호에서 N 번째 이후의 성분은

동일한 위상을 가지는 것을 알 수 있습니다.

 

따라서 정상적인 상황에서 CP 와 CP' 는 동일한 위상을 가져야합니다.

하지만 만약 반송파 주파수의 오프셋(CFO)이 생기면

 

매 데이터 마다 CFO 가 누적되어 한바퀴를 돌고나면

CFO * N 에 해당하는 Offset 이 발생합니다.

 

 

CFO*N 값은 두 신호중 하나를 아래로 뒤집어서 (conjugate)

곱한 결과의 위상 성분에 해당합니다.

따라서 해당 위상에

(-N pi) 를 나눠 주면 반송파의 Offset 을 알 수 있습니다.

 

그리고 CFO 를 보상해줄 때는

(-k pi) 를 곱한 값으로 회전(e의 승수를 곱해줌) 시켜주면 됩니다.

(K 는 해당 데이터 번호)

수식적인 유도

 

 

 

MATLAB 코드 1 (Van de Beek 논문 코드 구글링)

%% ML estimation of time and frequency offset in OFDM systems
% Algorithm authors : J.J. van de Beek, M. Sandell, and P.O. B훟rjesson
% Prepared by Hiren Gami

clc
clear all;
%% Parameter declaration according to paper
Nfft = 1024;
Ncp = 128;
Nsym = 6;
FreqOffset = 2.3;
SNRdb = 10;
theta = 256;
%% OFDM symbol generation
data = 2*randi([0 1],1, Nsym*Nfft)-1; % BPSK data

Tx = zeros(1,Nsym*(Nfft+Ncp));
OFDMsym = zeros(1,Nfft);
for sym = 1:Nsym
    OFDMsym = ifft(data(Nfft*(sym-1)+1:(Nfft*sym)),Nfft)*sqrt(Nfft); % data -> ifft sym (per sym)
    Tx((Nfft+Ncp)*(sym-1)+1:(Nfft+Ncp)*sym) = [OFDMsym(Nfft-Ncp+1:Nfft) OFDMsym]; % add CP
end

%% AWGN channel
snr = 10^(-SNRdb/10);
noise = sqrt(snr/2)*(randn(1,Nsym*(Nfft+Ncp))+1i*randn(1,Nsym*(Nfft+Ncp)));
Rx = exp(1i*2*pi*FreqOffset*(0:length(Tx)-1)./Nfft).*Tx + noise;
c = length(Tx);
%% ML estimation of timing and frequency offset
PHI_sum = zeros(1,Nsym*(Nfft+Ncp)-Nfft);
GM_sum = zeros(1,Nsym*(Nfft+Ncp)-Nfft);

for n = 1:Nsym*(Nfft+Ncp)-(Nfft+Ncp)
    PHI=0;GM=0;
    for m = n:n+Ncp-1
        PHI = PHI+ (Rx(m)*conj(Rx(m)) + Rx(m+Nfft)*conj(Rx(m+Nfft)));
        GM = GM+ Rx(m)*conj(Rx(m+Nfft));
    end
    PHI_sum(n) = abs(GM)- (snr/(snr+1))*PHI;
    GM_sum(n) = -angle(GM)/(2*pi);
end

%% Estimation results display

[~, timeOffset] = max(PHI_sum);
freqOffset = GM_sum(timeOffset);

subplot(2,1,1);plot(PHI_sum);title('Estimation of timing offset');grid on;
subplot(2,1,2);plot(GM_sum);title('Estimation of frequency offset');grid on;

[value,thetaEst]=findpeaks(PHI_sum,'minpeakdistance',Nfft);
disp(thetaEst);
disp(GM_sum(thetaEst));

MATLAB 코드 2 (PSS detection 에서 CFO estimation 을 추가한 코드)

 

만약 CFO estimation 을 안한 결과를 보고 싶으시면

상단 부근의 Initialize 영역에서

bCfoCpEstimation 을 false 로 하시면 됩니다 ;)

 

 

pssDetectWithCfoCpEstimation.zip

0.00MB

pNotPss 는 Pss 가 아닌 신호를 받을 때 correlation 이 threshold 보다 높은지 체크

pPss 는 루트 값이 맞는지, 시간 오프셋이 맞는지, correlation 이 threshold 보다 높은지 체크한 결과.

CP 를 이용해 CFO 를 조정하니 훨씬 좋은 결과가 나타남!!

clear all; 

close all;

%% initialize ==============================================================
symbolSz = 1024; % this is used as fft size

NumSubcarrier = 600;

cpLen = 72;

cellID = randi([0 503]);

IdN1 = ceil(cellID / 3); % SSS

IdN2 = rem(cellID, 3); % PSS

corrFrameSize = 72;

padLen = corrFrameSize/2;

startPoint = cpLen / 4;

freqOffset = 0.3;

%% Generate data

data = rand(NumSubcarrier, 1) + 1i*rand(NumSubcarrier, 1);
data2 = rand(NumSubcarrier, 1) + 1i*rand(NumSubcarrier, 1);

data = [data(1:NumSubcarrier/2); 0; data(NumSubcarrier/2+1 : end)];
data2 = [data2(1:NumSubcarrier/2); 0; data2(NumSubcarrier/2+1 : end)];

%% Generate PSS ===========================================================

zcRootI = [25, 29, 34];

pss = [];
for n = 1 : 3
    tmp = zadoffChuSeq(zcRootI(n), 63);
    pss = cat(2, pss, [zeros(5, 1);tmp(1:31); zeros(1, 1); tmp(33:end);zeros(5, 1)]);
end

pss_tx = pss(1 + IdN2 * 73 : (1 + IdN2) * 73);
pss_tx = pss_tx';

%% generate SSS

[d0, d5] = generate_sss(IdN1);
tmp = randi(2);
if (tmp ~= 0)
    sss = d0;
else
    sss = d1;
end

sss_tx = [zeros(1,5) sss zeros(1,5)];

%% ifft

% n point ifft
pss_tx_time = ifft(pss_tx, symbolSz);
sss_tx_time = ifft(sss_tx, symbolSz)';
data_tx_time = ifft(data, symbolSz);
data2_tx_time = ifft(data2, symbolSz);

% add CP
pss_tx_cp_time = [pss_tx_time(1:cpLen); pss_tx_time];
data_tx_cp_time = [data_tx_time(1:1+cpLen); data_tx_time];
data2_tx_cp_time = [data2_tx_time(1:1+cpLen); data2_tx_time];

E_N = length(pss_tx_cp_time);
Eavg_tx_time = sum(abs(pss_tx_cp_time) .^ 2)/length(E_N);

%pssTxCpCfoTime = pss_tx_cp_time;
pssTxCpCfoTime = generateCFO(pss_tx_cp_time, freqOffset, symbolSz);

pss_ifft = [];
for n = 0:2
    tmp = ifft(pss(1 + n * 73 : (1 + n) * 73)', symbolSz);
    pss_ifft =  [pss_ifft; tmp];
end

%% pss detection ==========================================================

threshhold = [0 0.001+0.0003*(1:6) inf];


pPss = zeros(length(threshhold), 1);
pNotPss = zeros(length(threshhold), 1);

iter_num = 1e+3;

snr_dB = 0;

Pd_error_offset_s = [];
Pe_error_offset_s = [];

for n = 1:length(threshhold)

    for iter = 1:iter_num
        %% pss sig

        
        corr_maxs = zeros(3,1);
        cors = zeros(corrFrameSize, 1);
        detect_offsets = zeros(3,1);
        for m = 0:2

            % padding for time offset
            rand_padding2 = zeros(padLen,1);

            %sss_tx_n_time = my_awgn(sss_tx_time, snr_dB);
            data_tx_cp_n_time = my_awgn(data_tx_cp_time, snr_dB);
            pssTxCpNoiseCFOTime = my_awgn(pssTxCpCfoTime, snr_dB);

            %pss_tx_n_pad_time = [sss_tx_n_time((1+symbol_sz-pad_len : symbol_sz)); pss_tx_n_time; data_tx_n_time(1 : pad_len)];            
            pss_tx_cp_n_pad_time = [pssTxCpNoiseCFOTime; data_tx_cp_n_time(1:corrFrameSize)];

            % detecting
            for l = 0 : corrFrameSize-1

                
                pss_tx_sample = pss_tx_cp_n_pad_time(l+1 + startPoint:l+symbolSz + startPoint);

                %pss_tx_n_fft = fft(pss_tx_sample);

                % FD
                %cor = abs(xcorr(pss_tx_n_fft(1:73), pss(1 + m * 73 : (1 + m) * 73)', 0));

                % TD
                tmp = xcorr(pss_tx_sample, pss_ifft(1 + m * symbolSz : (1 + m) * symbolSz), 0);
                cors(1+l) = (real(tmp)^2+imag(tmp)^2);
            end
            [corr_maxs(m + 1), detect_offsets(m + 1)] = max(cors);
        end

        [max_cor, detect_ID_2] = max(corr_maxs);
        detect_offset = detect_offsets(detect_ID_2);

        % offset == pad_len because time offset is pad_len
        if (threshhold(n) <= max_cor && ...
                detect_ID_2-1 == IdN2)
            if (detect_offset == cpLen - startPoint + 1)
                pPss(n) =  pPss(n) + 1;
            else
                Pd_error_offset_s = [Pd_error_offset_s detect_offset];
            end
        end

        %% rand sig

        
        corr_maxs = zeros(3,1);
        detect_offsets = zeros(3,1);
        for m = 0:2

            data_pad_time = [data_tx_cp_time; data2_tx_cp_time];

            data_tx_cp_n_time = my_awgn(data_pad_time, snr_dB);

            % padding for time offset
            padLen = symbolSz/4;

            rand_padding2 = zeros(padLen,1);

            

            % detecting
            for l = 0 : corrFrameSize-1

                rand_tx_sample = data_pad_time(l+1+startPoint:l+symbolSz+startPoint);

                % FD
                %cor = abs(xcorr(rand_n_fft(1:73), pss(1 + m * 73 : (1 + m) * 73)', 0));

                % TD
                tmp = xcorr(rand_tx_sample, pss_ifft(1 + m * symbolSz : (1 + m) * symbolSz), 0);
                cors(1+l) = (real(tmp)^2+imag(tmp)^2);
            end

            corr_maxs(m + 1) = max(cors);
        end


        [max_cor, detect_ID_2] = max(corr_maxs);
        %detect_offset = detect_offsets(detect_ID_2);

        % offset == pad_len because time offset is pad_len
        %detect_offset == pad_len && ...
        if (threshhold(n) <= max_cor)
            pNotPss(n) =  pNotPss(n) + 1;
        end
    end
end

pNotPss = pNotPss / iter_num;
pPss = pPss / iter_num;
plot(pNotPss, pPss, 'b--o');
title('ROC')
xlabel('p Not PSS')
ylabel('p PSS')


%% get energy
function [energy] = get_energy(signal)
    % Get the number of symbols
    N = length(signal);
    
    % Calculate Symbol Energy
    energy = sum(abs(signal) .^ 2)/N;
end


%% awgn noise with power
function [noise] = get_wgn_energy(energy, signal, SNR_dB)

    % Specify SNR in dB. Try setting various different value here and see how the result changes
    
    % Get the number of symbols
    N = length(signal);
    
    % Convert SNR (in dB) to SNR (in Linear)
    SNR_lin = 10 .^ (SNR_dB/10);
    
    % Calculate the Sigma (Standard Deviation) of AWGN
    awgnSigma = sqrt(energy/(2*SNR_lin));
     
    % Generate a sequence of noise with Normal Distribution and rescale it with the sigma
    noise = awgnSigma*(randn(1,N)+1i*randn(1,N));
    noise = noise';

end

function signalCFO = generateCFO(signal, cfo, fftSize)
    tmp =exp(1i*2*pi*cfo*(0:length(signal)-1)./fftSize);
    signalCFO = tmp'.*signal; 
end

function [timeOffset, freqOffset] = calculateCfoCp(signal, nFft, nCp, nSym, snrDb)

    PHI_sum = zeros(1,nSym*(nFft+nCp)-nFft);
    GM_sum = zeros(1,nSym*(nFft+nCp)-nFft);
    
    for n = theta:nSym*(nFft+nCp)-(nFft+nCp)
        PHI=0;GM=0;
        for m = n:n+nCp-1
            PHI = PHI+ (signal(m)*conj(signal(m)) + signal(m+nFft)*conj(signal(m+nFft)));
            GM = GM+ signal(m)*conj(signal(m+nFft));
        end
        PHI_sum(n) = abs(GM)- (snrDb/(snrDb+1))*PHI;
        GM_sum(n) = -angle(GM)/(2*pi);
    end
    
    [~, timeOffset] = max(PHI_sum);
    freqOffset = GM_sum(timeOffset);
end