Re: [Discuss-gnuradio] OFDM demodulation problem / example video

[Achilleas Anastasopoulos]

Jens,

I would like to elaborate on some of your comments.
I am using the following notation:
N fft size (2048)
G guard interval (504)
v0 fractional frequency offset
f0=v0 T absolute frequency offset
tau0 fractional time offset
t0=tau0 T absolute timne offset

> To sum up: 
>
> Frequency offset results in ICI. 
> 	This is "white noise" interference.

Correct: v0 will result in a time offset of N*v0
which will manifest itself as ICI.

There is MORE to it.
It ALSO introduces an unknown phase offest (constant in each subchannel):
Consider the i-th component of consequtive OFDM symbols.
They are effectively rotated by multiples of 2*pi*(N+G)*v0.
This is the phase accumulated over an OFDM symbol duration.
Although the linear phase increase will be compensated by D-QPSK, there 
will be a constant phase rotation equal to 2*pi*(N+G)*v0 in all
demodulated QPSK symbols.

BTW, If you make sure that you estimate v0 with an accuracy of 1e-6
then the time offset is on the order of 1e-3 (negligible) and the
frequency offset on the order of 1 degree (negligible).

BTW, because of the way you have coded your frequency correction (in 
your Matlab code) even perfect estimation of v0 will still result in the 
above phenomenon.
If you correct for v0 INDIVIDUALLY in each OFDM symbol (as you do in 
your code), then--assuming perfect estimation of v0--there will be two 
effects:
1) within each OFDM symbol there will be no phase variation since you 
have corrected it already
2) each consequtive sample in the i-th subchannel will be advanced by
2*pi*(N+G)*v0 radians. This is the phase accumulated over an OFDM symbol 
duration.
This is why it is important to derotate your entire frame at once.


Having said that, it seems that even if this is done, the problem you 
observe is still there!



> Time offset results in 
>
>         a) ISI, if offset results in taking values from other symbols.
>         This is also "white" interference.
>
>         b) phase shift, if offset results in taking values from the
>         cyclic prefix. This results in a "circle" structure, as the
>         phase shift is higher for subcarriers of higher frequency.
>
> The time offset effect b) is removed by differential demodulation, since
> the phase shift will be the same.

I agree with the above.
Also, the above assumes that the offset tau0 remains fixed for
each OFDM symbol. Otherwise, the comment made by Bob is valid.
However, I think the DAB system is designed to operate under this 
assumption, since they have a single pilot in the begining of the frame.

> The effect I am seeing is a constant phase shift for all subcarriers, 
> different for each OFDM symbol. Nothing that could be explained with a
> frequency offset or a time offset. 
>
> I am clueless.

I attach a piece of Matlab code that simulates both the Tx and Rx.
Everything mentioned above can be demonstrated in this controled experiment.

I still do not understand why the real-data application does not work.

Achilleas

function [f0 ,t0] = freq_time_est(ff,M,sync,wframe)



rec=zeros(M,2);
for t=1:M,
%t

out=zeros(1,length(ff));
for k=1:length(ff)
    f=ff(k);
    s=sync.*exp(j*2*pi*f*[1:length(sync)]);
    out(k)=mean(s.*conj(wframe(t:t+length(sync)-1)));
end
[value f0i]=max(out);

rec(t,:)=[value f0i];


end
[gvalue t0]=max(rec(:,1));
f0=ff(rec(t0,2));

clear all;
%close all;


Q=76;
%data=zeros(1537,Q-1);
data=floor(rand(1537,Q-1)*4);

qpsk=exp(j*2*pi*data/4);
qpsk(769,:)=zeros(1,Q-1);

%dqpsk modulation
dqpsk=zeros(1537,Q);
sync=floor(rand(1537,1)*4);
%dqpsk(:,1)=ones(1537,1);
dqpsk(:,1)=exp(j*2*pi*sync/4);
dqpsk(769,1)=0;
for k=2:Q,
    dqpsk(:,k)=dqpsk(:,k-1).*qpsk(:,k-1);
end

% initial phases of the oscilators
phase0=2*pi*rand(1537,1);
dqpsk=dqpsk.*(exp(j*phase0)*ones(1,Q));

% ofdm modulation
frame=zeros(504+2048,Q);
for l=1:Q
    a=dqpsk(:,l).';
    aa=[zeros(1,256) a zeros(1,255)];
    aaa=fftshift(aa);
    am=ifft(aaa); % OFDM symbol
    amp=[am(end-503:end) am].'; % with prefix
    amp=amp/sqrt(mean(abs(amp).^2)); % normalize power to 1
    frame(:,l)=amp;
end
msync=frame(:,1);  % modulated sync (1st ofdm symbol with prefix) for later use

D=1000;   
frame_h=[zeros(1,D) reshape(frame,1,Q*(504+2048)) zeros(1,D)]; % zeros in head 
and tail








%=======================================================
%channel
% ISI
% Le=300;
% h=randn(1,Le)+j*randn(1,Le);
% h=h./sqrt(mean(abs(h).^2));
% frame_h=conv(h,frame_h);

% frequency error
df=1.0e-5;
f=frame_h.*exp(j*2*pi*df*[1:length(frame_h)]+j*rand(1,1)*2*pi);


% timing error (uncomment this if you want to introduce timing error)
% M=11;   % max integer error [-M,M] samples
% ts0=0;
% while ts0==0
%    ts0=floor(rand(1,1)*M)-floor((M-1)/2);
% end
% ts0; % so the timing error is ts0 samples (so far) 
% 
% L=5;   % resolution of fractional errors
% tf0=0;
% while tf0==0
%    tf0=floor(rand(1,1)*L)-floor((L-1)/2);
% end
% tf0; % so the timing error is tt0/L samples (fractional)
% 
% tt0=ts0*L+tf0; % total timing error
% [ts0 tf0/L tt0/L]
% framei_h=resample(f,L,1);
% t0=D*L+1+tt0;
% rframe_h=framei_h(1:L:end); 
% rframe_h=rframe_h(1:D+Q*(504+2048)+D); % received signal



% %=======================================================
% %frequency/timing estimation (correlation with msync)
% % we do not use oversampled receiver (easy to extend)
% % coarse
% Lag=floor(M/2)+2;
% tframe=rframe_h(D+1-Lag:D+504+2048+Lag);  % you know the sync start is close 
to D+1 
% [f0 t0]=freq_time_est(-1e-2:1e-4:1e-2,2*Lag,msync.',tframe);
% [f0 t0]
% rframe_h=rframe_h.*exp(-j*2*pi*f0*[1:length(rframe_h)]);
% rframe_h=rframe_h(t0-Lag:end); % recenter so the new estimate of the sync 
begining is in D+1
% % fine
% Lag=5;
% tframe=rframe_h(D+1-Lag:D+504+2048+Lag);
% [f0 t0]=freq_time_est(-1e-4:1e-5:1e-4,2*Lag,msync.',tframe);
% [f0 t0]
% rframe_h=rframe_h.*exp(-j*2*pi*f0*[1:length(rframe_h)]);
% rframe_h=rframe_h(D+t0-Lag:D+t0-Lag+(504+2048)*Q-1); % sync starts here


% comment this line if you introduce timing errors
rframe_h=f(D+1:D+Q*(504+2048)); % received signal (correct timing)



rframe=reshape(rframe_h,504+2048,Q);

%demod
dqpsk_est=zeros(1537,Q);
for l=1:Q
    y = fftshift(fft(rframe(504+1:end,l),2048));   
    z = y(257:1793);
    z = z / sqrt(mean(abs(z).^2));
    dqpsk_est(:,l)=z;
end

qpsk_est=zeros(1537,Q-1);
data_est=zeros(1537,Q-1);
for l=1:Q-1
    dem=dqpsk_est(:,l+1).*conj(dqpsk_est(:,l));
    qpsk_est(:,l)=dem;
%     plot(real(dem), imag(dem), 'o');
%     axis([-3 3 -3 3]);
%     drawnow
%     pause
end

plot(qpsk_est(:,:).*exp(-j*0*pi/180),'o');grid on;
drawnow;

% calc effective snr
qq=reshape(qpsk_est,1,1537*(Q-1));
ii=find(real(qq)>0.5);
q=qq(ii);
n=q-mean(q);

ii=find(real(qq)<-0.5);
q=qq(ii);
n=[n q-mean(q)];

ii=find(imag(qq)>0.5);
q=qq(ii);
n=[n q-mean(q)];

ii=find(imag(qq)<-0.5);
q=qq(ii);
n=[n q-mean(q)];

snr=-10*log10(mean(abs(n).^2))