Fetching contributors…
Cannot retrieve contributors at this time
97 lines (86 sloc) 3.43 KB
% Multi-Frame Analysis based on derivative of Linear Interpolation (dLI-MFA)
% Input
% af : Array of structures with matrices containing sinusoidal
% parameters (as the output of sin_analysis.m).
% Each matrix is made of:
% The 1st line for the frequency of each sinusoid [Hz]
% The 2nd line for their amplitude (linear scale)
% The DC has to be included (in the first column)
% fs : [Hz] Signal's sampling frequency
% [extrap_dcny] : If true (default), extrapolate sinusoidal components at DC
% and up to Nyquist.
% If false, use the sinusoidal components as they are.
% Output
% E : The amplitude cepstral envelope
% Reference
% [1] G. Degottex, "A Time Regularization Technique for Discrete Spectral
% Envelopes Through Frequency Derivative", Signal Processing Letters, IEEE,
% 22(7):978-982, July 2015.
% Copyright (c) 2013 Foundation for Research and Technology-Hellas - Institute
% of Computer Science (FORTH-ICS)
% License
% This file is part of libphoni. libphoni is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser General
% Public License as published by the Free Software Foundation, either version 3
% of the License, or (at your option) any later version. libphoni is
% distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
% PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
% details.
% Author
% Gilles Degottex <>
function [E] = env_dli_mfa(af, fs, dftlen, extrap_dcny)
if nargin<4; extrap_dcny=true; end
% If asked, extrapolate components at DC and up to Nyquist
if extrap_dcny; af = env_extrap_sins_dcny(af, fs); end
% For each frame, compute the frequency derivative of the linear envelope
F = fs*(0:dftlen/2)/dftlen;
dA = ones(numel(af),dftlen/2);
for fi=1:numel(af)
A = interp1(af(fi).sins(1,:), log(af(fi).sins(2,:)), F, 'linear', 'extrap');
dA(fi,:) = diff(A); % derivative approximation
% Smooth the frequency derivative across time
win = hamming(size(dA,1));
win = win./sum(win);
dAw = dA.*repmat(win, 1, size(dA,2));
mA = sum(dAw,1);
% Retrieve the envelope
E = cumsum(mA);
E = [E(1), E];
% Align the envelope on the harmonics of the central frame
ci = floor((numel(af)-1)/2)+1;
ek = interp1(F', E, af(ci).sins(1,:));
ak = log(af(ci).sins(2,:));
idx = find(af(ci).sins(1,:)>0 & af(ci).sins(1,:)<4000);
E = E + mean(ak(idx)) - mean(ek(idx));
E = exp(E);
% Plot the final solution
if 0
hold off;
plot(F(1:end-1), dA, 'r');
hold on;
plot(F(1:end-1), mA, 'b');
xlim([0 fs/2]);
xlabel('Frequency [Hz]');
hold off;
plot(F, mag2db(E), 'b');
hold on;
plot(af(ci).sins(1,:), mag2db(af(ci).sins(2,:)), 'xk');
xlim([0 fs/2]);
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]');
% keyboard