% function ofdm_tx(N_data)
function yv = ofdm_tx(N_data)

j = sqrt(-1);
WORK_BANDPASS = 1;
do_evaluate = 0;

% Params
fa = 48000;
fc = fa/4;
F1 = 750;
F2 = 2250;
F3 = 3000;

N_Carrier = 226;
M_qam = 64;
N_lp = 301;
TSym = 21.3330E-3
TGrd = 5.333E-3
TFrm = TSym+TGrd
NGrd = round(TGrd*fa)
NSym = round(TSym*fa)
NFrm = round(TFrm*fa)
N_fft = NSym

bin_c = fix(N_fft*fc/fa)+1;
bin_f1 = fix(N_fft*(fc+F1)/fa)+1;
bin_f2 = fix(N_fft*(fc+F2)/fa)+1;
bin_f3 = fix(N_fft*(fc+F3)/fa)+1;

% Source
qlut = qamtable(M_qam);
s = round(M_qam*rand(N_data,1)+0.5);

% Band pass for transmission
f_lp = N_Carrier*fa/N_fft
hbp = cos(2*pi*(0:N_lp-1)*fc/fa).*sinc((-(N_lp-1)/2:(N_lp-1)/2)*f_lp/fa);
hbp = hbp/sum(hbp.^2).*Hanning(N_lp)';
hlp = sinc((-(N_lp-1)/2:(N_lp-1)/2)*f_lp/fa);
hlp = hlp/sum(hlp.^2).*Hanning(N_lp)';

t = 1;
if (WORK_BANDPASS == 1)
    bin_range = N_fft/4-N_Carrier/2:N_fft/4+N_Carrier/2;
else
    bin_range = 1:N_Carrier;
end;

N_bins = length(bin_range);
for k=1:N_data/N_bins,

    % DC
    X = zeros(N_fft,1);
    X(1) = 0;

    for m=1:N_bins,
        % Data
        X(bin_range(m)) = qlut(s(t));
        t = t + 1;
    end;

    % Additional carrier
    if (WORK_BANDPASS == 1)
        X(bin_c) = 1.0;
        X(bin_f1) = 1.0 + j*1.0;
        X(bin_f2) = 1.0 - j*1.0;
        X(bin_f3) = -1.0 + j*1.0;
    end;

    % Make conj. symmetric in F => real in T
    X(N_fft:-1:N_fft/2+1) =  conj(X(2:N_fft/2+1));

    % Modulate
    x = ifft(X);
    xgrd = x(NSym-NGrd+1:NSym);
    xt = [xgrd' x']';
    xv((k-1)*(NFrm)+1:k*(NFrm)) = xt;
    
    if (do_evaluate == 1)
        xf = fft(xt(NGrd+1:NSym+NGrd), N_fft);
        plot(xf, '+'); grid;
        pause;
    end;    

end;

k_fft = 0.3333*sqrt(N_fft);

% Pulse shape
if (WORK_BANDPASS == 1)
    yv = real(filter(hbp,1,k_fft*xv));
   % yv = k_fft*xv;
else
    yc = cos(2*pi*(0:length(xv)-1)*fc/fa);
    yF1 = 1.0/N_fft*cos(2*pi*(0:length(xv)-1)*F1/fa);
    yF2 = 1.0/N_fft*cos(2*pi*(0:length(xv)-1)*F2/fa);
    yF3 = 1.0/N_fft*cos(2*pi*(0:length(xv)-1)*F3/fa);
    yv = yc.*(0.025 + filter(hlp,1,k_fft*xv) + yF1 + yF2 + yF3);
end;

% Output
close all;

figure;
freqz(hlp);

wavwrite(yv,fa,16,'yv.wav')