ita_result2audio.m 5.03 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
function varargout = ita_result2audio(varargin)
%ITA_RESULT2AUDIO - converts itaResult to itaAudio by interpolation and resampling
%  This function converts an itaResult into an itaAudio by interpolation
%  and resampling
%
%  Syntax:
%   audioObjOut = ita_result2audio(resultObjIn,samplingRate,fftDegree,options)
%   Options (default):
%    'no_filter' (false) :       
%
%  Example:
%   audioObjOut = ita_result2audio(resultObjIn,44100)
%
%   Reference page in Help browser 
%        <a href="matlab:doc ita_result2audio">doc ita_result2audio</a>

% <ITA-Toolbox>
% This file is part of the ITA-Toolbox. Some rights reserved. 
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder. 
% </ITA-Toolbox>


% Author: MMT -- Email: mmt@akustik.rwth-aachen.de
% Created:  19-Jan-2010 

%% Get Function String
thisFuncStr  = [upper(mfilename) ':'];     %Use to show warnings or infos in this functions

%% Initialization and Input Parsing
% if a target samplingrate and fftdegree is specified
if nargin < 3
    if nargin < 2
        varargin = [varargin, {44100}];
    end
    varargin = [varargin, {10}];
end

sArgs        = struct('pos1_input','itaResult','pos2_targetSamplingRate','numeric','pos3_fftDegree','numeric','no_filter',false);
[input,targetSamplingRate,fftDegree,sArgs] = ita_parse_arguments(sArgs,varargin); 

%% 'input' is an audioObj and is given back
% first do the interpolation to a minFFTDegree
minFFTDegree = nextpow2(input.nSamples);
if isTime(input)
    oldAbscissa  = input.timeVector;
    samplingRate = 1./mean(diff(oldAbscissa));
    minFFTDegree = max(minFFTDegree,log2(ceil(max(oldAbscissa)*samplingRate/2)*2));
    newAbscissa  = (0:2^minFFTDegree-1)./samplingRate;
else
    oldAbscissa  = input.freqVector;
    samplingRate = ceil(2*max(oldAbscissa)); % round to 1Hz precision
    minFFTDegree = max(minFFTDegree,log2(ceil(samplingRate/2/mean(diff(oldAbscissa)))));
    newAbscissa  = linspace(0,samplingRate/2,(2^minFFTDegree)+1);
end
oldData = input.data;
interpData = zeros(numel(newAbscissa),size(oldData,2));
for i=1:input.nChannels
    if all(isreal(oldData(:,i)))
        interpData(:,i) = interp1(oldAbscissa,real(oldData(:,i)),newAbscissa,'spline','extrap');
    else % let's try abs and unwrapped phase (or groupdelay)
        interpData(:,i) = interp1(oldAbscissa,abs(oldData(:,i)),newAbscissa,'spline','extrap').*...
            exp(1i.*interp1(oldAbscissa,unwrap(angle(oldData(:,i))),newAbscissa,'spline','extrap'));
%         interpData(:,i) = interp1(oldAbscissa,real(oldData(:,i)),newAbscissa,'spline','extrap')+...
%             1i.*interp1(oldAbscissa,imag(oldData(:,i)),newAbscissa,'spline','extrap');
    end
end

% save some values for postprocessing
processingLimits = [min(oldAbscissa) max(oldAbscissa)];
% avoid the zero
if isFreq(input)
    processingLimits(1) = find(oldAbscissa~=0,1,'first');
end

% if an fftdegree was specified that is higher than the
% minFFTDegree, do another interpolation
if fftDegree > minFFTDegree
    ita_verbose_info([thisFuncStr 'fftDegree was specified, doing second interpolation'],1);
    oldAbscissa = newAbscissa;
    if isTime(input)
        newNSamples = max(2^fftDegree,ceil(max(oldAbscissa)*samplingRate/2)*2);
        newAbscissa = (0:newNSamples-1)./samplingRate;
    else
        newAbscissa = linspace(0,samplingRate/2,2^(fftDegree-1)+1);
    end
    oldData = interpData;
    interpData = zeros(numel(newAbscissa),size(oldData,2));
    for i=1:size(oldData,2)
        if all(isreal(oldData(:,i)))
            interpData(:,i) = interp1(oldAbscissa,real(oldData(:,i)),newAbscissa,'spline','extrap');
        else % let's try abs and unwrapped phase (or groupdelay)
            interpData(:,i) = interp1(oldAbscissa,abs(oldData(:,i)),newAbscissa,'spline',1e-6).*...
                exp(1i.*interp1(oldAbscissa,unwrap(angle(oldData(:,i))),newAbscissa,'spline',1e-6));
        end
    end
end

if ~isreal(interpData(1,:))
    interpData(1,:) = real(interpData(1,:));
end
if ~isreal(interpData(end,:))
    interpData(end,:) = real(interpData(end,:));
end
sObj = saveobj(input);
sObj = rmfield(sObj,[{'classname','classrevision'},itaResult.propertiesSaved]);
input = itaAudio(sObj);
input.signalType = 'energy';
input.samplingRate = samplingRate;
input.data = interpData;

% some postprocessing to cancel out effects of extrapolation
if isTime(input)
113 114 115 116 117 118
    if processingLimits(1) == 0
        %only use fadeOut part of timeWindow
        input = ita_time_window(input,[0.8*processingLimits(2) 0.95*processingLimits(2)],'time');
    else
        input = ita_time_window(input,[1.1*processingLimits(1) 0.95*processingLimits(1) 0.8*processingLimits(2) 0.95*processingLimits(2)],'time');
    end
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
elseif ~sArgs.no_filter
    input = ita_filter_bandpass(input,'lower',processingLimits(1),'upper',processingLimits(2),'order',20,'zerophase');
end

if targetSamplingRate ~= samplingRate
    input = ita_resample(input,targetSamplingRate);
    input = ita_interpolate_spk(input,fftDegree);
end

%% Add history line
input = ita_metainfo_add_historyline(input,mfilename,varargin);

%% Set Output
varargout(1) = {input}; 

%end function
end