%demo           Script to analyse the speech signal shown in section IV.D of
%               the paper and to compute the harmonic component using
%               different envelope modeling techniques.
%
%This has been tested and works using MATLAB 5.3, Signal Processing Toolbox 4.2
%and Optimization Toolbox 2.0. It should work with earlier versions of these
%granted that you replace aryule by lpc (in the main script) and fminunc by
%fmin (in the function OLC).
%
%_______________________________________________________________________
%For more information about this mfile please take a look at the paper
%"Estimation of the spectral envelope of voiced sounds using a penalized
%likelihood approach" by Marine Oudot, Olivier Cappé and Eric Moulines.

% Author: Olivier Cappé, June 2000
% ENST Dpt. TSI / LTCI (CNRS URA 820), Paris

%%% 1. Analysis parameters
% Sampling frequency (Hz)
smpfreq = 8000;
% Half the frame length
frmext = 120; % 15 ms
% Hop size
hopsze = 40; % 5 ms
% Data taper
taper = hanning(2*frmext+1);
st = sum(taper);
st2 = sum(taper.^2);
% Optimization options for OLC
options = optimset('LargeScale', 'off', 'GradObj', 'on', 'TolFun', 1e-3, ...
  'TolX', 1e-3, 'MaxFunEvals', 250, 'LineSearchType', 'cubicpoly', ...
  'Display', 'off', 'DerivativeCheck', 'off');
% Frequency axis for plotting spectra
wtr = linspace(0, pi, 1024);

%%% 2. Read signal and pitch track
% Load input signal (signal stored in IEEE little endian short format)
fid = fopen('signal.raw', 'rb',  'ieee-le');
sig = fread(fid, 'short');
fclose(fid);
% Load pitch table (time in s, frequency in Hz) and convert to sample and
% radian frequency
load signal.dat;
pitch = 2*pi*signal(:,2)/smpfreq;
% Position of the analysis frames (use the same as for the pitch)
frmloc =  1 + round(smpfreq*signal(:,1));
frmnb = length(frmloc);

%%% 3. Declare arrays
% Synthetic signals
synth_ar = zeros(size(sig));
synth_olc = zeros(size(sig));
% Normalization
olanorm = zeros(size(sig));

%%% 4. Main analysis loop (i is the frame index)
for i = 1:frmnb
  % Only do something when the frame is voiced
  if (~isnan(pitch(i)))
    %%% 4.a. Compute harmonic amplitudes and apprent noise level
    % Signal frame
    ind = (frmloc(i)-frmext:frmloc(i)+frmext)';
    s = sig(ind);
    % Compute number of harmonics
    hrmnb = floor(pi*0.96/pitch(i));
    % Compute the amplitudes of the harmonics using tapered Fourier transform
    REG = exp(j*pitch(i)*ind*(1:hrmnb));
    theta = REG'*(taper.*s)*(2/st);
    a = abs(theta).^2;
    % Estimate noise psd
    cn = robnest(s-real(REG*theta), pitch(i), 16);
    % Apparent noise level
    n = (4*st2/(st^2))*exp(cepsval(pitch(i)*(1:hrmnb), cn));

    %%% 4.b. Estimation of the envelope using AR
    [a_ar, e_ar] = aryule(s.*taper, 12);
    % Energy correction for windowing
    e_ar = e_ar*((2*frmext+1)/st2)*2/hrmnb;
    
    %%% 4.c. Estimation using OLC or WLS
    % Preliminary estimate using discrete power cepstrum
    c_ls = LS(a, pitch(i), 40);
    % Estimate the enveloppe using OLC
    c_olc = OLC(c_ls, a, pitch(i)*(1:hrmnb), n, options);
    % Alternatively use WLS by uncommenting this line (and commenting the
    % previous one)
    %c_olc = WLS(a, pitch(i)*(1:hrmnb), 40, n);
    
    %%% 4.d. Plots
    % Compute the apparent nois epsd estimate
    nhat = (4*st2/(st^2))*exp(cepsval(wtr, cn));
    % Compute enveloppe estimated using AR
    ehat_ar = 20*log10(abs(freqz(sqrt(e_ar), a_ar, wtr)));
    % Compute enveloppe estimated using OLC
    ehat_olc = 10*cepsval(wtr, c_olc)/log(10);
    clf;
    S = fft(s.*taper, 256)*(2/st);
    plot(2*pi*(0:128)/256, 20*log10(abs(S(1:129)))); % Periodogram
    hold on;
    plot(wtr, 10*log10(nhat), ':'); % estimated noise level
    plot(wtr, ehat_ar, 'k'); % AR envelope
    plot(wtr, ehat_olc, '--r'); % OLC or WLS envelope
    xlabel('radian frequency');
    ylabel('dB');
    title(['Frame ' int2str(i)]);
    drawnow;
    
    %%% 4.e. Signal reconstruction using OLA
    % Compute the harmonic magnitudes from the envelope
    mamp_ar = abs(freqz(sqrt(e_ar), a_ar, pitch(i)*(1:hrmnb)));
    mamp_olc = exp(cepsval(pitch(i)*(1:hrmnb), c_olc)/2);
    % Time domain synthetic signal using modeled amplitudes and original phases
    phase = angle(theta);
    s_ar = real(REG*(mamp_ar'.*exp(j*phase)));
    s_olc = real(REG*(mamp_olc'.*exp(j*phase)));
    % OLA
    synth_ar(ind) = synth_ar(ind) + s_ar.*taper;
    synth_olc(ind) = synth_olc(ind) + s_olc.*taper;
    % Keep track of normalization
    olanorm(ind) = olanorm(ind) + taper;
  end
end
%%% 5. Normalize the signals
ind = find(olanorm > 0);
synth_ar(ind) = synth_ar(ind)./olanorm(ind);
synth_olc(ind) = synth_olc(ind)./olanorm(ind);
