function fdet_wav(filename, foff, p)
%
% Example: fdet_wav('drm-15435-if.wav', -174.65, [700 850 1])
% Example: ofdm_tx(40, 1)
%          fdet_wav('xv.wav', 0, [700 800 1])
%       or fdet_wav('yv.wav', 0, [700 800 1])

close all;
fs = 48000;
fc = fs/4;
N = 1024;
Ng = 256;

omega_min = p(1)/fs;
omega_max = p(2)/fs;
omega_step = p(3)/fs;

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

[len ndim] = size(xwav);

if (ndim > 1)
    fprintf('Mode: Baseband\n');
    x = xwav(:, 1) + j*xwav(:, 2);
else
    fprintf('Mode: Passband\n');
    x = xwav.*exp(-j*(2*pi*fc/fs*(0:len-1)'));
end

fprintf('Perform frequency correction\n');
x = x.*exp(j*(2*pi*foff/fs*(0:len-1)'));

% Remove cyclic prefix
fprintf('Perform CP calculation\n');
x0 = x(1:2*(N+Ng));
x1 = x(N+Ng+1:3*(N+Ng));

c = xcorr(x0, x1);
[v, ii] = max(abs(c));
cp_off = mod(ii, Ng+N)

% Detect CFO
fprintf('Perform frequency detection\n');
[omega, err, dmod, omega_coarse, omega_est, Xm] = fdet(x, N, Ng, cp_off, omega_min, omega_step, omega_max);

f_coarse = omega_coarse*fs
f_est = omega_est*fs

fprintf('Receive Symbols\n');
K = fix((length(x)-Ng)/(N+Ng));
for k=10:30
    xp = x((k-1)*(N+Ng)+cp_off+1:k*(N+Ng)+cp_off);
    Xp(:,k) = fft(xp(Ng+1:N+Ng), N);
end;

close all;

plot(abs(c)); grid;
figure;

bin = fix(N*750/48000)+1
plot(real(Xp(bin, :)), imag(Xp(bin, :)), '+'); grid;

figure;
plot(0:N-1, abs(Xm)); grid;

figure;
subplot(5,1,1)
plot(1:lge(omega), real(dmod)/N, 1:lge(omega), imag(dmod)/N, 1:lge(omega), abs(dmod)/N); grid; legend('Re', 'Im', 'abs()');

subplot(5,1,2)
plot(1:lge(omega), err*fs, '-'); grid; legend('error');

subplot(5,1,3)
plot(1:lge(omega), omega*fs); grid; legend('f');

subplot(5,1,4)
plot(1:lge(omega), angle(dmod), '-', 1:lge(omega), diff(angle([0 dmod])), 'r-'); grid; legend('phi', 'dphi');

subplot(5,1,5)
plot(1:lge(omega), log10(abs(dmod) ./ (1E-6+abs(real(dmod)))), '-'); grid; legend('diff');