clear all;
close all;

% dab.dat contains data sampled from the USRP 2 MHz (decimation factor 32)

fid = fopen('dab_ml.dat', 'r');
%skip 50000 samples - bug in gnuradio?
d = fread(fid, 50000, 'float32');

%read real
d_re = fread(fid, 1000000, 'float32',4);
fclose(fid);

fid = fopen('dab_ml.dat', 'r');
%skip 50001 samples - bug in gnuradio?
d = fread(fid, 50001, 'float32');

%read imag
d_im = fread(fid, 1000000, 'float32', 4);
fclose(fid);

e = d_re + j*d_im;
clear d*

% not all data need, file too big
chunk = e(1:700000);

% Generate phase reference symbol, frequency domain
p = prs;

% resampling
chunk_resamp = resample(chunk, 2048, 2000);

% this code is used to look for the synch frame manually
if 1==0

glob_max = 0;
for k = 3.6*10^5:length(chunk_resamp)
    data = chunk_resamp(k:k+2048);
    f = fftshift(fft(data,2048));

    xc = abs(xcorr(p,f));

    [m offset] = max(xc);
    if m > 10^7
        disp('found synch frame');
        plot(xc);
        axis([0 4000 0 2*10^7]);
        drawnow
        pause
        k
    end
    if m > glob_max
        offset_max = offset
        glob_max = m
        glob_max_k = k
        plot(xc);
        axis([0 4000 0 2*10^7]);
        drawnow
        %pause
    end    
    drawnow
end

end

% decode one frame (75 symbols + synch symbol)
for o = 0:75

% start of frame, found manually
start_resamp = 168988;
    
%get the data for one OFDM symbol (2048 samples), skip guard intervall (504
%samples)
data_resamp = chunk_resamp(start_resamp+(504+2048)*o:start_resamp+2048-1+(504+2048)*o);

% compensate frequency offset - done manually
data2_resamp = data_resamp.*exp(j*linspace(0,2*pi*(0.385-8),length(data_resamp)).');

% fft to decode data
d_resamp = fftshift(fft(data2_resamp,2048));

% only 1536 subcarriers + DC carrier to display
data = d_resamp(256:1792);

% first OFDM symbol is used for synchronization - no data.
if o == 0
    plot(abs(d_resamp).^2);    
    drawnow
    pause
end

% display decoded data symbols, including DC carrier
if o ~= 0
    dem = data ./ data_l;

    plot(real(dem), imag(dem), 'o');   
    axis([-3 3 -3 3]);

    drawnow
   
    pause
end

% DAB uses pi/4 DPSK, last frame is reference
data_l = data;

    
end