function [sym_out, H, W] = ofdm_rx(filename, foff, flip_spec, SNR, with_cir, show_H)
%
% [sym_out, H, W] = ofdm_rx(filename, foff, flip_spec, SNR, with_cir, show_H)
% Example: ofdm_rx('drm-15435-if.wav', -174.65, 0, 400, 0, 1)
% Example: ofdm_rx('xv.wav', 0, 0, 400, 0, 1)
%       or ofdm_rx('yv.wav', 0, 0, 400, 0, 1)

close all;

drm_mode = 'B';
drm_bw = '10k';

[ofdm_drm_params, ofdm_spec_occ] = drm_params(drm_mode, drm_bw);

pre_mix = 0;
fs = 12000;
N = 256;
Ng = 64;
[carrier_pilot, phi_pilot, mag_pilot] = getRefFreq(drm_mode);
[carrier_time, phi_time, mag_time] = getRefTime(drm_mode);

h_cir = [0 0 .1 -.2 .5 .2 -.1 0 0 ];

% Generate signal
[xwav fs_ nbits] = wavread(filename);

[len ndim] = size(xwav);

if (ndim == 2)
    fprintf('Mode: Baseband\n');
    x = xwav(:, 1) + j*xwav(:, 2);
    x = x.*exp(-j*(2*pi*(foff/fs_)*(0:len-1)'));
elseif (ndim == 1)
    fprintf('Mode: Passband\n');
	x = xwav - mean(xwav);
    x = xwav.*exp(-j*(2*pi*(0.25-foff/fs_)*(0:len-1)'));
    hlp = FIRCalcLowpass(12000/fs_, 15);
    x = filter(hlp, 1, x);
	
	% Downsample
	x = x(1:4:length(x));
else
    error('Cannot read file!');
end

if with_cir
    x = filter(h_cir, 1, x);
end

x = x - mean(x);
if flip_spec
	x = conj(x);
end
x = x./max(abs(x));
xout = [real(x) imag(x)];
wavwrite(xout,fs,nbits,'ofdm_rx.wav');

% Modulate on carrier
xup = upsample(x, 4);
hlp = FIRCalcLowpass(0.25, 129);
xup = filter(hlp, 1, xup);

y = exp(j*2*pi.*(0:length(xup)-1)'*0.25) .* xup;

wavwrite(y,4*fs,nbits,'ofdm_rx_passband.wav');

x = awgn(x, SNR, 'measured');
K = fix((length(x))/(N+Ng))-1;

% Timing coarse
fprintf('Coarse timing acquisition\n');
cp_off = 0;
cp_off_ = [];
Z = 0;
for k=1:20
	xp = x((k-1)*(N+Ng)+1:k*(N+Ng));
	[cp_off, Nd, valid, Z] = timing_coarse(xp, N, Ng, Ng, Z);
	if (valid)
		cp_off_ = [cp_off_ cp_off];
	end
end;
cp_off = round(median(cp_off_))
Nd = Nd

% Integer CFO coarse
fprintf('Coarse integer CFO acquisition\n');
fdet_coarse_Z = 0;
Z_cfa_pre = 0;
df_acq_coarse_int_ = [];
valid_latency_count = 20;
for k=1:200
    xp = x((k-1)*(N+Ng)+cp_off+1:k*(N+Ng)+cp_off);
    [d_bin, valid, fdet_coarse_Z] = cfo_coarse_int(xp(Ng+1:N+Ng), N, carrier_pilot, 11, fdet_coarse_Z);
    if valid
        df_acq_coarse_int = -d_bin(1)/N*fs;
        if valid_latency_count > 0
            valid_latency_count = valid_latency_count - 1;
        else
            break;
        end
    end
end

fprintf('Coarse fractional CFO acquisition\n');
% CFO coarse
df_acq_coarse_fract_ = [];
for k=1:200
    xp = x((k-1)*(N+Ng)+cp_off+1:k*(N+Ng)+cp_off);
    [domega_f, Z_cfa_pre] = cfo_coarse_fract(xp, N, Ng, Ng, Nd, Z_cfa_pre);
    df_acq_coarse_fract_ = [df_acq_coarse_fract_ domega_f*fs];
end;
df_acq_c = df_acq_coarse_int + mean(df_acq_coarse_fract_(10:200))
df_acq_c_ = df_acq_coarse_int + df_acq_coarse_fract_;

if pre_mix
    x1 = x.*exp(j*(2*pi*df_acq_c/fs.*(0:length(x)-1)'));
    df_acq_fine = 0;
else
    x1 = x;
    df_acq_fine = df_acq_c;
end

% CFO Fine
fprintf('Fine fractional CFO acquisition\n');
Xlast = [];
df_acq_fine_err = 0;
df_acq_fine_err_ = [];
df_acq_fine_ = [];
kn = fs/(2*pi*(1+Ng/N))/N;
x_ = [];
Hp_ = [];
Z_ = [];
phi = 0.0;
bin_track = 1;
sym_pilot = exp(j*2*pi*phi_pilot'/1024);

bin_pilot = c2i(N, carrier_pilot);
for k=1:K,
    xp = x1((k-1)*(N+Ng)+cp_off+1:k*(N+Ng)+cp_off);
    xp = xp.*exp(j*(2*pi*df_acq_fine/fs*(0:N+Ng-1)' + phi));
    phi = phi + 2*pi*df_acq_fine/fs*(N+Ng);
    xp = xp(Ng+1:N+Ng);
    x_ = [x_' xp']';
    X = fft(xp);
    if (~isempty(Xlast))
        Z = Xlast .* conj(X);
        df_acq_fine_err = kn*angle(Z(bin_pilot));
        df_acq_fine_err_ = [df_acq_fine_err_ df_acq_fine_err];
        phi_p = angle(X(bin_pilot));
        mag_p = abs(X(bin_pilot));
        df_acq_fine = df_acq_fine + 0.1*mean(df_acq_fine_err(bin_track));
        df_acq_fine_ = [df_acq_fine_ df_acq_fine];
        Hp = X(bin_pilot)./conj(sym_pilot);
        Hp_ = [Hp_ Hp];
        Z_ = [Z_ X];
    end
    Xlast = X;
end

figure;
plot(0:length(cp_off_)-1, cp_off_, 0:length(cp_off_)-1, cp_off*ones(size(cp_off_)), 'r-'); grid; title('cp_{off}')

df_acq_fine_ = df_acq_fine_';

figure;
subplot(3, 1, 1)
plot(0:length(df_acq_c_)-1, df_acq_c_); grid; title('df_{coarse} (Acquisition)')

subplot(3, 1, 2)
plot(0:length(df_acq_fine_)-1, df_acq_fine_); grid; title('df_{fine} (Acquisition)')

subplot(3, 1, 3)
plot(0:length(df_acq_fine_err_)-1, df_acq_fine_err_); grid; title('Error(df_{fine}) (Acquisition)')

figure;
Z_pilot = Z_(bin_pilot, 200:length(Z_))';
plot(real(Z_pilot), imag(Z_pilot), '+'); grid; title('SymPilot_{Freq}'); 
maxZ = max(max(abs(Z_pilot)));
if (maxZ < 1)
    axis([-1 1 -1 1]);
else
    axis([-maxZ maxZ -maxZ maxZ]);
end

fprintf('Frame timing acquisition\n');
%% Detect start of DRM frame
bins = c2i(N, carrier_time);
kc_max = -1000;
n_max = 0;
skip_syms = 15*4;
for n=0:15-1
    frame_timing_Z = 0;
    k = skip_syms+n;
    for m=1:2
        [kc valid frame_timing_Z] = timing_frame(Z_(bins, k)', bins, 15, frame_timing_Z);
        k = k + 15;
    end
    if valid
        if kc > kc_max
            n_max = n;
            kc_max = kc;
        end
    end
end
frame_start = n_max

figure;
start = skip_syms+frame_start;
Z_timing = Z_(bins, start:15:length(Z_))';
plot(real(Z_timing), imag(Z_timing), '+'); grid; title('SymPilot_{Time}'); 
maxZ = max(max(abs(Z_timing)));
if (maxZ < 1)
    axis([-1 1 -1 1]);
else
    axis([-maxZ maxZ -maxZ maxZ]);
end

fprintf('Channel estimation and symbol reception\n');

[W_syms, W_pilots] = calcWiener(drm_mode, drm_bw);

figure;
% Detect scattered gain pilots
symbols_per_frame = ofdm_drm_params.nspf;
symbols_to_delay = ofdm_drm_params.y;
frames_per_window = 2*ofdm_drm_params.y;
H_gain_cell = cell (frames_per_window, 1);
H_gain_carriers_cell = cell (frames_per_window, 1);
H_raw = zeros(length(ofdm_spec_occ.kmin:ofdm_spec_occ.kmax),1);
carriers = ofdm_spec_occ.kmin:ofdm_spec_occ.kmax;
carrier_indexes = c2i(N, carriers);
fac_cell_list = ...
{
    [],
    [],
    [13, 25, 43, 55, 67],
    [15, 27, 45, 57, 69],
    [17, 29, 47, 59, 71],
    [19, 31, 49, 61, 73],
    [9, 21, 33, 51, 63, 75],
    [11, 23, 35, 53, 65, 77],
    [13, 25, 37, 55, 67, 79],
    [15, 27, 39, 57, 69, 81],
    [17, 29, 41, 59, 71, 83],
    [19, 31, 43, 61, 73],
    [21, 33, 45, 63, 75],
    [23, 35, 47, 65, 77],
    []
};
s = 0;
m = 0;

H_syms = [];
H_pilots = [];
sym_out = [];
kk = 1;
H_ = [];
df_track_ = [];
MIN_ABS_H = 1e-10;

for k=start:length(Z_)
    
    %shifted symbol index
    n = rem(m+symbols_to_delay-1, ofdm_drm_params.y);

	[carrier_gain, phi_gain, mag_gain] = getRefGain(ofdm_drm_params, ofdm_spec_occ, s);
    carrier_gain_bins = c2i(N, carrier_gain);
	H_  = transpose(Z_(carrier_gain_bins, k))./(mag_gain.*exp(j*2*pi*phi_gain/1024));

    H_raw(carrier_gain_bins) = (H_);
    H_gain_cell{m+1} = (H_);
    H_gain_carriers_cell{m+1} = carrier_gain;
    
    str = sprintf('Frame %d, s=%d, m=%d, n=%d', k, s, m, n);
    if kk > 6
		ii = 1+mod(m+symbols_to_delay-1+(0:frames_per_window-1), frames_per_window);
		H_gain = [H_gain_cell{ii}];
		H_syms = H_gain*W_syms{n+1};
		H_pilots = (H_gain)*W_pilots{n+1};
		ZZ = Z_(carrier_indexes, k-symbols_to_delay).';
		hhh = H_gain_cell{mod(m-1, frames_per_window)+1};
		df_track = 12000/(2*pi*N)*angle( (H_pilots) * hhh' + MIN_ABS_H);
		df_track_ = [df_track_ df_track];
		
		% check for too small values
		rms_H = sqrt( abs( H_syms*H_syms' )/length(H_syms) );
		MIN_ABS_H = max( rms_H/20, 1e-10 );	%quick and dirty fixed settings
		H_too_small_index = find( abs(H_syms)<MIN_ABS_H );
		H_syms( H_too_small_index ) = MIN_ABS_H;

		sym_eq = ZZ ./ H_syms;
		sym_out(:, kk) = sym_eq;
		if 0
			subplot(2, 1, 1)
			plot(1:length(hhh), abs(hhh), 'b-', 1:length(H_pilots), abs(H_pilots), 'r-'); grid; 
			subplot(2, 1, 2)
			plot(1:length(hhh), angle(hhh), 'b-', 1:length(H_pilots), angle(H_pilots), 'r-'); grid; 
			title(str);
			pause(0.01);
		end
		if 1
        if show_H
			H_raw(1:2:N) = H_raw(2:2:N);
			subplot(3, 1, 1)
			plot(carriers, abs(H_raw(carrier_indexes)), 'b-', carriers, abs(H_syms), 'r-'); grid; title('abs(H_{syms})'); axis([ofdm_spec_occ.kmin ofdm_spec_occ.kmax 0 15+0*max(abs(H_syms))]); title(str);
			subplot(3, 1, 2)
			plot(carriers, angle(H_raw(carrier_indexes)), 'b-', carriers, angle(H_syms), 'r-'); grid; title('phi(H_{syms})'); axis([ofdm_spec_occ.kmin ofdm_spec_occ.kmax -4 4]); title(str);
			subplot(3, 1, 3)
			plot(carriers, abs(sym_eq), 'b-'); grid; title('abs(Sym_{Eq})'); axis([ofdm_spec_occ.kmin ofdm_spec_occ.kmax 0 2]); title(str);
            hold on;
			for b=c2i(N, carrier_pilot)
	            x = [b b];
		        y = [0 2];
			    plot(x-1,y, 'r-');
			end
            hold off;

            pause(0.01);
		else
            fac_c = fac_cell_list{mod(s-symbols_to_delay, symbols_per_frame)+1};
            if (~isempty(fac_c))
				fac_i = fac_c - ofdm_spec_occ.kmin + 1;
				hold on;
		        plot(sym_eq(fac_i), '+', 'MarkerSize', 4);
			    pause(0.01);
				grid; title('FAC'); axis([-1 1 -1 1]);
				hold off;
			end
		end
		end
    end
    kk = kk + 1;
    m = mod(m+1, frames_per_window);
    s = mod(s + 1, symbols_per_frame);

end

figure;
plot(df_track_); grid; title('df_{track}');

function Y = spec_order(N, X)
    Y = [X(N/2+1:N)' X(1:N/2)']';