SebastianEwert@28: function [f_audio_out,timepositions_afterDegr] = degradationUnit_adaptiveEqualizer(f_audio, samplingFreq, timepositions_beforeDegr, parameter)
SebastianEwert@28: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SebastianEwert@28: % Name: degradationUnit_adaptiveEqualizer
SebastianEwert@28: % Date of Revision: 2013-01-23
SebastianEwert@28: % Programmer: Sebastian Ewert
SebastianEwert@28: %
SebastianEwert@28: % Description:
SebastianEwert@28: % - designs a filter such that the mean spectrum of f_audio becomes similar
SebastianEwert@28: % to a given mean spectrum
SebastianEwert@28: % - mean spectra are specified using a decibel scale, i.e. if x is a magnitude
SebastianEwert@28: % spectrum, then apply x -> 20*log10(x)
SebastianEwert@28: % - there are four ways to specify the destination mean spectrum: 1. by
SebastianEwert@28: % loading example files provided with the toolbox, 2. by using specific
SebastianEwert@28: % noise "color" profiles, 3. by providing the destination mean spectrum
SebastianEwert@28: % using the parameter destMagFreqResp, 4. by providing audio data from
SebastianEwert@28: % which the destination mean spectrum is computed.
SebastianEwert@28: % - if mean spectra are computed, this is done by computing a magnitude
SebastianEwert@28: % spectrogram, deriving the corresponding values in dB, and then averaging
SebastianEwert@28: % over all time frames.
SebastianEwert@28: % - the same is done for f_audio and the two mean spectral vectors are
SebastianEwert@28: % compared
SebastianEwert@28: % - Then a filter is designed such that f_audio's mean spectral vector
SebastianEwert@28: % becomes similar to destMagFreqResp
SebastianEwert@28: % - filtering is done using a linear-phase FIR filter 
SebastianEwert@28: % - this unit normalizes the output
SebastianEwert@28: %
SebastianEwert@28: % Notes:
SebastianEwert@28: % - note that the mean spectrum is often a good approximation of what is
SebastianEwert@28: % sometimes referred to as the mean spectral shape (however, this term is not
SebastianEwert@28: % defined for general sound recordings). In this sense, the function
SebastianEwert@28: % modifies the mean spectral shape to the desired one.
SebastianEwert@28: %
SebastianEwert@28: % Input:
SebastianEwert@28: %   f_audio      - audio signal \in [-1,1]^{NxC} with C being the number of
SebastianEwert@28: %                  channels
SebastianEwert@28: %   samplingFreq - sampling frequency of f_audio
SebastianEwert@28: %   timepositions_beforeDegr - some degradations delay the input signal. If
SebastianEwert@28: %                             some points in time are given via this
SebastianEwert@28: %                             parameter, timepositions_afterDegr will
SebastianEwert@28: %                             return the corresponding positions in the
SebastianEwert@28: %                             output. Set to [] if unavailable. Set f_audio
SebastianEwert@28: %                             and samplingFreq to [] to compute only
SebastianEwert@28: %                             timepositions_afterDegr.
SebastianEwert@28: %
SebastianEwert@28: % Input (optional): parameter
SebastianEwert@28: %   .loadInternalMagFreqResp=1  - loads one of the destMagFreqResp provided
SebastianEwert@28: %                                 by the toolbox
SebastianEwert@28: %   .internalMagFreqResp='Beatles_NorwegianWood'
SebastianEwert@28: %   .computeMagFreqRespFromAudio - computes destMagFreqResp from given
SebastianEwert@28: %                                  audio data
SebastianEwert@28: %   .computeMagFreqRespFromAudio_audioData
SebastianEwert@28: %                               - audio data for .computeMagFreqRespFromAudio
SebastianEwert@28: %   .computeMagFreqRespFromAudio_sf - sampl freq for .computeMagFreqRespFromAudio
SebastianEwert@28: %   .destMagFreqResp = []       - in db. See above.
SebastianEwert@28: %   .destMagFreqResp_freqs = [] - must have same length as destMagFreqResp.
SebastianEwert@28: %                                 In Hertz
SebastianEwert@28: %   .fftLength = 2 ^ nextpow2(0.02 * samplingFreq); - fft length to
SebastianEwert@28: %                                 calculate spectrogram of f_audio.
SebastianEwert@28: %
SebastianEwert@28: % Output:
SebastianEwert@28: %   f_audio_out   - audio signal \in [-1,1]^{NxC} with C being the number
SebastianEwert@28: %                   of channels
SebastianEwert@28: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SebastianEwert@28: 
SebastianEwert@28: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SebastianEwert@28: % Audio Degradation Toolbox
SebastianEwert@28: %
SebastianEwert@28: % Centre for Digital Music, Queen Mary University of London.
SebastianEwert@28: % This file copyright 2013 Sebastian Ewert, Matthias Mauch and QMUL.
SebastianEwert@28: %
SebastianEwert@28: % This program is free software; you can redistribute it and/or
SebastianEwert@28: % modify it under the terms of the GNU General Public License as
SebastianEwert@28: % published by the Free Software Foundation; either version 2 of the
SebastianEwert@28: % License, or (at your option) any later version.  See the file
SebastianEwert@28: % COPYING included with this distribution for more information.
SebastianEwert@28: %
SebastianEwert@28: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SebastianEwert@28: 
SebastianEwert@28: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SebastianEwert@28: % Check parameters
SebastianEwert@28: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SebastianEwert@28: if nargin<4
SebastianEwert@28:     parameter=[];
SebastianEwert@28: end
SebastianEwert@28: if nargin<3
SebastianEwert@28:     timepositions_beforeDegr=[];
SebastianEwert@28: end
SebastianEwert@28: if nargin<2
SebastianEwert@28:     error('Please specify input data');
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: if isfield(parameter,'loadInternalMagFreqResp')==0
SebastianEwert@28:     parameter.loadInternalMagFreqResp = 1;
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'loadNoiseColorPreset')==0
SebastianEwert@28:     parameter.loadNoiseColorPreset = 0;
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'computeMagFreqRespFromAudio')==0
SebastianEwert@28:     parameter.computeMagFreqRespFromAudio = 0;
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: if isfield(parameter,'internalMagFreqResp')==0
SebastianEwert@28:     parameter.internalMagFreqResp = 'Beethoven_Appasionata_Rwc';
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'noiseColorPreset')==0
SebastianEwert@28:     parameter.noiseColorPreset = 'pink';
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'computeMagFreqRespFromAudio_audioData')==0
SebastianEwert@28:     parameter.computeMagFreqRespFromAudio_audioData = [];
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'computeMagFreqRespFromAudio_sf')==0
SebastianEwert@28:     parameter.computeMagFreqRespFromAudio_sf = [];
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'destMagFreqResp')==0
SebastianEwert@28:     parameter.destMagFreqResp = [];
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'destMagFreqResp_freqs')==0
SebastianEwert@28:     parameter.destMagFreqResp_freqs = [];
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'fftLength')==0
SebastianEwert@28:     parameter.fftLength = 2 ^ nextpow2(0.02 * samplingFreq);
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'filterOrder')==0
SebastianEwert@28:     parameter.filterOrder = round(parameter.fftLength/2);
SebastianEwert@28: end
SebastianEwert@28: if isfield(parameter,'visualizations')==0
SebastianEwert@28:     parameter.visualizations = 0;
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SebastianEwert@28: % Main program
SebastianEwert@28: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SebastianEwert@28: if isempty(f_audio)
SebastianEwert@28:     % we design a linear-phase filter. Such filters have the property that
SebastianEwert@28:     % the group delay for every frequency is always equal to
SebastianEwert@28:     % parameter.filterOrder/2. Therefore, although we built a signal
SebastianEwert@28:     % adaptive filter, we can adjust for that delay. That means that the
SebastianEwert@28:     % timepositions_beforeDegr don't require any adjustment.
SebastianEwert@28:     timepositions_afterDegr = timepositions_beforeDegr;
SebastianEwert@28:     return;
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: % load/compute destMagFreqResp
SebastianEwert@28: [destMagFreqResp,destMagFreqResp_freqs] = internal_getDestMagFreqResp(parameter);
SebastianEwert@28: 
SebastianEwert@28: % compute mean spectral vector for f_audio
SebastianEwert@28: [meanMagSpec,meanMagSpec_freqs] = internal_computeMeanSpectralVector(f_audio,samplingFreq,parameter.fftLength);
SebastianEwert@28: meanMagSpec = internal_standardizeMagFreqResp(meanMagSpec);
SebastianEwert@28: 
SebastianEwert@28: % compute magnitude response for the filter to be designed
SebastianEwert@28: destMagFreqResp_org = destMagFreqResp;
SebastianEwert@28: destMagFreqResp_freqs_org = destMagFreqResp_freqs;
SebastianEwert@28: if ~((length(destMagFreqResp_freqs)==length(meanMagSpec_freqs)) && (all(meanMagSpec_freqs(:) == destMagFreqResp_freqs(:))))
SebastianEwert@28:     % in this case we interpolate the frequency response using a
SebastianEwert@28:     % spline interpolation
SebastianEwert@28:     destMagFreqResp = spline(destMagFreqResp,destMagFreqResp_freqs,meanMagSpec_freqs);
SebastianEwert@28: end
SebastianEwert@28: filter_magFreqResp = destMagFreqResp(:) - meanMagSpec(:);
SebastianEwert@28: 
SebastianEwert@28: % design filter (fir2 is linear phase)
SebastianEwert@28: filter_magFreqResp_linear = 10 .^ (filter_magFreqResp/20);
SebastianEwert@28: b = fir2(parameter.filterOrder,meanMagSpec_freqs/(samplingFreq/2),filter_magFreqResp_linear);
SebastianEwert@28: 
SebastianEwert@28: % apply filter
SebastianEwert@28: parameterApplyImpulseResponse.loadInternalIR = 0;
SebastianEwert@28: parameterApplyImpulseResponse.impulseResponse = b;
SebastianEwert@28: parameterApplyImpulseResponse.impulseResponseSampFreq = samplingFreq;
SebastianEwert@28: parameterApplyImpulseResponse.normalizeOutputAudio = 1;
SebastianEwert@28: parameterApplyImpulseResponse.averageGroupDelayOfFilter = round(parameter.filterOrder/2);
SebastianEwert@28: [f_audio_out,timepositions_afterDegr] = degradationUnit_applyImpulseResponse(f_audio, samplingFreq, timepositions_beforeDegr, parameterApplyImpulseResponse);
SebastianEwert@28: 
SebastianEwert@28: if parameter.visualizations
SebastianEwert@28:     fvtool(b,1);
SebastianEwert@28:     
SebastianEwert@28:     [meanMagSpecOut,meanMagSpecOut_freqs] = internal_computeMeanSpectralVector(f_audio_out,samplingFreq,parameter.fftLength);
SebastianEwert@28:     meanMagSpecOut = internal_standardizeMagFreqResp(meanMagSpecOut);
SebastianEwert@28:     
SebastianEwert@28:     figure;
SebastianEwert@28:     plot(destMagFreqResp_freqs_org,destMagFreqResp_org,'y');
SebastianEwert@28:     hold on;
SebastianEwert@28:     plot(meanMagSpecOut_freqs,meanMagSpecOut,'k');
SebastianEwert@28:     title('Comparison: destMagFreqResp(y) and mean spectral vector of output(k)')
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: 
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: function [f_meanmagspec_db,freqs] = internal_computeMeanSpectralVector(f_audio,fs,fftLength)
SebastianEwert@28: 
SebastianEwert@28: f_audio = mean(f_audio,2);
SebastianEwert@28: 
SebastianEwert@28: [f_spec,freqs,time] = spectrogram(f_audio,hanning(fftLength),fftLength/2,fftLength,fs);
SebastianEwert@28: 
SebastianEwert@28: f_magspec_db = 20 * log10(abs(f_spec));
SebastianEwert@28: 
SebastianEwert@28: f_magspec_db(:,isinf(sum(abs(f_magspec_db),1))) = []; % ignore columns with -inf/inf entries
SebastianEwert@28: f_magspec_db(:,isnan(sum(abs(f_magspec_db),1))) = [];
SebastianEwert@28: 
SebastianEwert@28: f_meanmagspec_db = mean(f_magspec_db,2);
SebastianEwert@28: 
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: function magFreqResp = internal_standardizeMagFreqResp(magFreqResp)
SebastianEwert@28: 
SebastianEwert@28: temp = magFreqResp(~isinf(magFreqResp));
SebastianEwert@28: temp = temp(~isnan(magFreqResp));
SebastianEwert@28: maxRobust = max(temp);
SebastianEwert@28: 
SebastianEwert@28: magFreqResp = magFreqResp - maxRobust;
SebastianEwert@28: 
SebastianEwert@28: magFreqResp(magFreqResp > 0) = 0;  % remaining positive inf
SebastianEwert@28: magFreqResp(magFreqResp < -80) = -80;  % remaining positive inf
SebastianEwert@28: 
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: function [destMagFreqResp,destMagFreqResp_freqs,fftLength] = internal_getDestMagFreqResp(parameter)
SebastianEwert@28: if parameter.loadInternalMagFreqResp
SebastianEwert@28:     % load example included in toolbox
SebastianEwert@28:     
SebastianEwert@28:     fullFilenameMfile = mfilename('fullpath');
SebastianEwert@28:     [pathstr,name,ext] = fileparts(fullFilenameMfile);
SebastianEwert@28:     dirRootIRs = fullfile(pathstr,'../degradationData');
SebastianEwert@28:     
SebastianEwert@28:     names_internal = {'Beatles_NorwegianWood','Beethoven_Appasionata_Rwc'};
SebastianEwert@28:     indexInternal = find(strcmpi(names_internal,parameter.internalMagFreqResp), 1);
SebastianEwert@28:     if isempty(indexInternal)
SebastianEwert@28:         error('Please specify a valid internal name')
SebastianEwert@28:     end
SebastianEwert@28:     
SebastianEwert@28:     switch indexInternal
SebastianEwert@28:         case 1
SebastianEwert@28:             file = fullfile(dirRootIRs,'SpecEnvelopes/Beatles_NorwegianWood.mat');
SebastianEwert@28:         case 2
SebastianEwert@28:             file = fullfile(dirRootIRs,'SpecEnvelopes/Beethoven_Appasionata_Rwc.mat');
SebastianEwert@28:     end
SebastianEwert@28:     load(file, 'destMagFreqResp', 'destMagFreqResp_freqs');
SebastianEwert@28:     
SebastianEwert@28: elseif parameter.loadNoiseColorPreset
SebastianEwert@28:     switch(lower( parameter.noiseColorPreset))
SebastianEwert@28:         case 'white'
SebastianEwert@28:             freqExponent = 0;
SebastianEwert@28:         case 'pink'
SebastianEwert@28:             freqExponent = 0.5;
SebastianEwert@28:         case 'brown'
SebastianEwert@28:             freqExponent = 1;
SebastianEwert@28:         case 'blue'
SebastianEwert@28:             freqExponent = -0.5;
SebastianEwert@28:         case 'violet'
SebastianEwert@28:             freqExponent = -1;
SebastianEwert@28:     end
SebastianEwert@28:     
SebastianEwert@28:     lengthMagResp = parameter.fftLength/2+1;
SebastianEwert@28:     destMagFreqResp_freqs = linspace(0,samplingFreq/2,lengthMagResp);
SebastianEwert@28:     magResp = 1./destMagFreqResp_freqs.^freqExponent;
SebastianEwert@28:     magResp(1) = 1;
SebastianEwert@28:     destMagFreqResp = 20 * log10(magResp);
SebastianEwert@28:     
SebastianEwert@28: elseif parameter.computeMagFreqRespFromAudio
SebastianEwert@28:     % compute destMagFreqResp as mean spectral vector from given audio data
SebastianEwert@28:     [destMagFreqResp,destMagFreqResp_freqs] = internal_computeMeanSpectralVector(...
SebastianEwert@28:         parameter.computeMagFreqRespFromAudio_audioData,parameter.computeMagFreqRespFromAudio_sf,parameter.fftLength);
SebastianEwert@28: else
SebastianEwert@28:     destMagFreqResp = parameter.destMagFreqResp;
SebastianEwert@28:     destMagFreqResp_freqs = parameter.destMagFreqResp_freqs;
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: % standardize destMagFreqResp
SebastianEwert@28: destMagFreqResp = internal_standardizeMagFreqResp(destMagFreqResp);
SebastianEwert@28: end
SebastianEwert@28: 
SebastianEwert@28: 
SebastianEwert@28: 
SebastianEwert@28: 
SebastianEwert@28: 
SebastianEwert@28: