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

close all;

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

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

fs = 12000;
N = 256;
Ng = 64;

[carrier_pilot, phi_pilot, mag_pilot] = getRefFreq(drm_mode);
[carrier_time, phi_time, mag_time] = getRefTime(drm_mode);

% Read wav file
[fsize fs4] = wavread(filename, 'size')
[xwav, fs4, nbits] = wavread(filename, 0);
wav_numCh = fsize(2);
wav_size = fsize(1);

isModePassband = 0;
if (wav_numCh == 1)
    fprintf('Mode: Passband\n');
	omega_xwav = (0.25-foff/fs4);
	isModePassband = 1;
elseif (wav_numCh == 2)
    fprintf('Mode: Baseband\n');
	omega_xwav = (foff/fs4);
else
	error('Cannot read file');
end

buf_recv = buffer();
buf_Z = buffer();

fileChunkSize = N;
fileCunkCounter = 0;
phi_xwav = 0;
filter_h = FIRCalcLowpass(0.25, 15);
filter_z = [];

fileRemain = wav_size;

% Timing coarse stuff
timingCoarse = struct('Z', 0, 'sync_counter', 20, 'sync', 0, 'cp_off', 0, 'cp_off_', [], 'Nd', 0);

% CFO coarse integer stuff
cfoCoarseInteger = struct('Z', 0, 'sync_counter', 20, 'sync', 0, 'df', 0);

% CFO coarse fractional stuff
cfoCoarseFractional = struct('Z', 0, 'sync_counter', 200, 'sync', 0, 'err', 0, 'err_', [], 'df', 0, 'df_', []);

% CFO Fine stuff
cfoFine = struct('Z', 0, 'sync_counter', 200, 'sync', 0, 'err', 0, 'err_', [], 'df', 0, 'df_', [], 'Xlast', [], 'phi', 0, 'bin_track', 1);

% Frame timing stuff
frameTiming = struct('Z', 0, 'sync_counter', 20, 'sync', 0, 'kc_max', -1000, 'n_max', 0);

% Equalization stuff
MIN_ABS_H = 1e-10;
equalizer = struct('Z', 0, 'sync_counter', 20, 'sync', 0, 's', 0, 'm', 0, 'kk', 1, 'eps_abs_h', MIN_ABS_H);

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

H_syms = [];
H_pilots = [];
H_ = [];
df_track_ = [];

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_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],
    []
};

while(fileRemain > fileChunkSize)
	
	n0 = fileCunkCounter*fileChunkSize+1;
	n1 = (fileCunkCounter+1)*fileChunkSize;
	[xwav] = wavread(filename, [n0 n1]);
	
	if (~isModePassband)
	    xwav = xwav(:, 1) + j*xwav(:, 2);
	end
	xwav = xwav - mean(xwav);

	x = xwav.*exp(-j*(2*pi*omega_xwav*(0:fileChunkSize-1)' + phi_xwav)); 
	phi_xwav = mod(phi_xwav + 2*pi*omega_xwav*fileChunkSize, 2*pi);
	
	if flip_spec
		x = conj(x);
	end

	% add noise
	x = awgn(x, SNR, 'measured');
	
	if (isModePassband)
		[x, filter_z] = filter(filter_h, 1, x, filter_z);
		buf_recv.write(x(1:4:length(x)));
	else
		buf_recv.write(x);
	end
	
	% --------------------------------------------------------------
	% Timing coarse
	if (~timingCoarse.sync)
		x = buf_recv.read(N+Ng);
		[cp_off, Nd, valid, timingCoarse.Z] = timing_coarse(x, N, Ng, Ng, timingCoarse.Z);
		if (valid)
			timingCoarse.cp_off_ = [timingCoarse.cp_off_ cp_off];
			if (timingCoarse.sync_counter > 0)
				timingCoarse.sync_counter = timingCoarse.sync_counter - 1;
			else
				timingCoarse.sync = 1;
			end
		end;
		if (timingCoarse.sync)
			timingCoarse.cp_off = round(median(timingCoarse.cp_off_));
			timingCoarse.Nd = Nd;
			Nd = timingCoarse.Nd
			cp_off = timingCoarse.cp_off
			
			% Reset buffer to Cp offset
			buf_recv.setReadIndex(timingCoarse.cp_off);
		end
	end
	
	% Integer CFO coarse
	if (timingCoarse.sync) && (~cfoCoarseInteger.sync)
		x = buf_recv.read(N+Ng);
		if (~isempty(x))
			[d_bin, valid, cfoCoarseInteger.Z] = cfo_coarse_int(x(Ng+1:N+Ng), N, carrier_pilot, 11, cfoCoarseInteger.Z);
			if valid
				cfoCoarseInteger.df = -d_bin(1)/N*fs;
				if (cfoCoarseInteger.sync_counter > 0)
					cfoCoarseInteger.sync_counter = cfoCoarseInteger.sync_counter - 1;
				else
					cfoCoarseInteger.sync = 1;
				end
			end
		end
	end
	
	% Fractional CFO coarse
	if (cfoCoarseInteger.sync) && (~cfoCoarseFractional.sync)
		x = buf_recv.read(N+Ng);
		if (~isempty(x))
			[err, cfoCoarseFractional.Z] = cfo_coarse_fract(x, N, Ng, Ng, Nd, cfoCoarseFractional.Z);
			cfoCoarseFractional.err_ = [cfoCoarseFractional.err_ err*fs];

			if (cfoCoarseFractional.sync_counter > 0)
				cfoCoarseFractional.sync_counter = cfoCoarseFractional.sync_counter - 1;
			else
				cfoCoarseFractional.sync = 1;
				cfoCoarseFractional.df = cfoCoarseInteger.df + mean(cfoCoarseFractional.err_);
				cfoCoarseFractional.df_ = cfoCoarseInteger.df + cfoCoarseFractional.err_;
				cfoFine.df = cfoCoarseFractional.df;
			end
		end
	end

	% CFO Fine
	if (cfoCoarseFractional.sync)
		kn = fs/(2*pi*(1+Ng/N))/N;
		sym_pilot = exp(j*2*pi*phi_pilot'/1024);
		bin_pilot = c2i(N, carrier_pilot);
		
		x = buf_recv.read(N+Ng);
		if (~isempty(x))
			x = x.*exp(j*(2*pi*cfoFine.df/fs*(0:N+Ng-1)' + cfoFine.phi));
			cfoFine.phi = cfoFine.phi + 2*pi*cfoFine.df/fs*(N+Ng);
			X = fft(x(Ng+1:N+Ng));
			if (~isempty(cfoFine.Xlast))
				Z = cfoFine.Xlast .* conj(X);
				cfoFine.err = kn*angle(Z(bin_pilot));
				cfoFine.err_ = [cfoFine.err_ cfoFine.err];
				cfoFine.df = cfoFine.df + 0.1*mean(cfoFine.err(cfoFine.bin_track));
				cfoFine.df_ = [cfoFine.df_ cfoFine.df];
				buf_Z.write(X);
			end
			cfoFine.Xlast = X;
		end

		bins = c2i(N, carrier_time);
		skip_syms = ofdm_drm_params.nspf*4;
		% Frame timing acquisition
		if (~frameTiming.sync) && (buf_Z.len() >= (N*(ofdm_drm_params.nspf+skip_syms)))
			buf_Z.setReadIndex(N*skip_syms);
			for n=0:ofdm_drm_params.nspf-1
				frameTiming.Z = 0;
				k = skip_syms+n;
				Z = buf_Z.read(N);
				if ~isempty(Z)
					for m=1:2
						[kc valid frameTiming.Z] = timing_frame(Z(bins)', bins, ofdm_drm_params.nspf, frameTiming.Z);
						k = k + ofdm_drm_params.nspf;
					end
					if valid
						if kc > frameTiming.kc_max
							frameTiming.n_max = n;
							frameTiming.kc_max = kc;
						end
					end
				end
			end
			if ~isempty(Z)
				buf_Z.setReadIndex(N*(frameTiming.n_max));
				frameTiming.sync = 1;
				frame_start = frameTiming.n_max
			end
		end

		if (frameTiming.sync)
			
			Z = transpose(buf_Z.read(N));
			if ~isempty(Z)
				
				equalizer.n = rem(equalizer.m + symbols_to_delay-1, ofdm_drm_params.y);

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

				H_raw(carrier_gain_bins) = (H_);
				H_gain_cell{equalizer.m+1} = (H_);

				str = sprintf('Frame %d, s=%d, m=%d, n=%d', equalizer.kk, equalizer.s, equalizer.m, equalizer.n);

				if equalizer.kk > frames_per_window

					ii = 1+mod(equalizer.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{equalizer.n+1};
					H_pilots = (H_gain)*W_pilots{equalizer.n+1};

					% Get the correctly delayed symbols
					Z = buf_Z.readAt(N, -N*symbols_to_delay);
					sym = Z(carrier_indexes).';

					% Get the correctly delayed raw channel estimate
					Hr = H_gain_cell{mod(equalizer.m-1, frames_per_window)+1};

					% get frequency error
					% correlate wiener filtered channel with raw estimate
					df_track = fs/(2*pi*N)*angle(H_pilots * Hr' + equalizer.eps_abs_h);
					df_track_ = [df_track_ df_track];

					% check for too small values
					rms_H = sqrt( abs( H_syms*H_syms' )/length(H_syms) );
					equalizer.eps_abs_h = max( rms_H/20, MIN_ABS_H );	%quick and dirty fixed settings
					H_too_small_index = find( abs(H_syms)<equalizer.eps_abs_h );
					H_syms( H_too_small_index ) = equalizer.eps_abs_h;

					% equalize symbol
					sym_eq = sym ./ H_syms;
					if 0
						subplot(2, 1, 1)
						plot(1:length(Hr), abs(Hr), 'b-', 1:length(H_pilots), abs(H_pilots), 'r-'); grid; 
						subplot(2, 1, 2)
						plot(1:length(Hr), angle(Hr), '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(equalizer.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
				equalizer.kk = equalizer.kk + 1;
				equalizer.m = mod(equalizer.m + 1, frames_per_window);
				equalizer.s = mod(equalizer.s + 1, symbols_per_frame);
			end
		end
	end
	% --------------------------------------------------------------
	fileRemain = fileRemain - fileChunkSize;
	fileCunkCounter = fileCunkCounter + 1;
end

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

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

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

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

if 0
xout = [real(x) imag(x)];
wavwrite(xout,fs,nbits,'rx_x.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,fs4,nbits,'rx_y.wav');

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)']';