matlab/ofdm/ofdm_rx.m

function [x] = ofdm_rx(filename, foff, knoise_dB, with_cir)
%
% Example: ofdm_rx('drm-15435-if.wav', -174.65, -400, 0)
% Example: ofdm_tx(40, 1)
%          ofdm_rx('xv.wav', 0, -400, 0)
%       or ofdm_rx('yv.wav', 0, -400, 0)

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.*exp(-j*(2*pi*(0.25-foff/fs_)*(0:len-1)'));
    hlp = FIRCalcLowpass(12000/fs_, 33);
    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);
x = x./max(abs(x));
xout = [real(x) imag(x)];
wavwrite(xout,fs,nbits,'ofdm_rx.wav');

x = x +  sqrt(0.5)*10^(knoise_dB/20)*(randn(size(x)) + j*randn(size(x)));


K = fix((length(x))/(N+Ng))-1;

% Timing coarse
cp_off = 0;
cp_off_ = [];
if 1
    fprintf('Coarse timing estimation\n');
    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
end

fdet_coarse_Z = 0;
Z_cfa_pre = 0;
df_i_ = [];
valid_latency_count = 20;
if 1
fprintf('Coarse integer CFO estimation\n');
% Integer CFO coarse
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_i = -d_bin(1)/N*fs;
        if valid_latency_count > 0
            valid_latency_count = valid_latency_count - 1;
        else
            break;
        end
    end
end;
end

fprintf('Coarse fractional CFO estimation\n');
% CFO coarse
df_f_ = [];
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_f_ = [df_f_ domega_f*fs];
end;
df_coarse = df_i + mean(df_f_(10:200))
df_coarse_ = df_i + df_f_;

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

% CFO Fine
fprintf('Fine CFO estimation\n');
fprintf('Receive Symbols\n');
Xlast = [];
df_fine = 0;
df_fine_ = [];
df_track_ = [];
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_track/fs*(0:N+Ng-1)' + phi));
    phi = phi + 2*pi*df_track/fs*(N+Ng);
    xp = xp(Ng+1:N+Ng);
    x_ = [x_' xp']';
    X = fft(xp);
    if (~isempty(Xlast))
        Z = Xlast .* conj(X);
        df_fine = kn*angle(Z(bin_pilot));
        df_fine_ = [df_fine_ df_fine];
        phi_p = angle(X(bin_pilot));
        mag_p = abs(X(bin_pilot));
        df_track = df_track + 0.1*mean(df_fine(bin_track));
        df_track_ = [df_track_ df_track];
        Hp = X(bin_pilot)./conj(sym_pilot);
        Hp_ = [Hp_ Hp];
        Z_ = [Z_ X];
    end
    Xlast = X;
end

df_track = mean(df_track_(10:200))
mean_df_fine = mean(df_fine_, 2)
mu_Hp =  mean(abs(Hp_'))'
sigma_Hp = var(abs(Hp_'))'
SNR_Hp = -10*log10(sigma_Hp)'

G = (conj(Hp)*1)./(abs(Hp).^2 * 1 + sigma_Hp)
Hp_1 = 1./Hp

if ~isempty(cp_off_)
    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}')
end

df_track_ = df_track_';

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

subplot(3, 1, 2)
plot(0:length(df_track_)-1, df_track_); grid; title('df_{track}')

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

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

%% Detect start of DRM frame
bins = c2i(N, carrier_time);
kc_max = -1000;
n_max = 0;
for n=0:15-1
    frame_timing_Z = 0;
    k = 200+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 = 200+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

figure;
% Detect scattered gain pilots
sym_count = 0;
start = 200+frame_start;
H = zeros(N, 1);
sym_H = [];
for k=start:length(Z_)
    [carrier_gain, phi_gain, mag_gain] = getRefGain(ofdm_spec_occ, ofdm_drm_params, sym_count);
    carrier_gain_bins = c2i(N, carrier_gain);
    Hp = Z_(carrier_gain_bins, k)./(mag_gain.*exp(-j*2*pi*phi_gain/1024))';
    H(carrier_gain_bins) = Hp;
    sym_count_str = sprintf('%d', sym_count);
    if 0
        plot(abs(Hp)); grid; title(sym_count_str);
        pause(0.1);
    end

    if (sym_count < 14)
        sym_count = sym_count+1;
    else
        H(1:2:N) = H(2:2:length(H));
        if 1
%            G = ones(size(H));
%            G = 1./(H+1e-6);
            G = (conj(H)*1)./(abs(H).^2 * 1 + sigma_Hp(1));
            subplot(4, 1, 1)
            plot(abs(spec_order(N, H))); grid; title('abs(H)');
            subplot(4, 1, 2)
            plot(angle(spec_order(N, H))); grid; title('phi(H)');
            subplot(4, 1, 3)
            plot(abs(spec_order(N, G))); grid; title('abs(G)');
            subplot(4, 1, 4)
            plot(angle(spec_order(N, G))); grid; title('phi(G)');
            pause(0.01);

            bin = c2i(N, 13);
            sym_H = [sym_H Z_(bin, k).*G(bin)];
        end
        H_gather = zeros(N, 1);
        sym_count = 0;
    end;
end

figure;
plot(real(sym_H'), imag(sym_H'), '+'); grid; title('Sym_{H}'); axis([-1.5 1.5 -1.5 1.5]);

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