function fdet_wav(filename, foff, p)
%
% Example: fdet_wav('drm-15435-if.wav', [900 1000 10])
% Example: ofdm_tx(40, 1)
%          fdet_wav('xv.wav', [700 800 1])
%       or fdet_wav('yv.wav', [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

x = x.*exp(j*(2*pi*foff/fs*(0:len-1)'));
[omega, err, dmod, omega_coarse, omega_est] = fdet(x, N, Ng, omega_min, omega_step, omega_max);

f_coarse = omega_coarse*fs
f_est = omega_est*fs

close all;
subplot(4,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(4,1,2)
plot(1:lge(omega), err*fs, '-'); grid; legend('error');

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

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