matthiasm@0: function [f_audio_out,timepositions_afterDegr] = degradationUnit_applyImpulseResponse(f_audio, samplingFreq, timepositions_beforeDegr, parameter) matthiasm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matthiasm@0: % Name: degradation_applyImpulseResponse matthiasm@0: % Date of Revision: 2013-04 matthiasm@0: % Programmer: Sebastian Ewert matthiasm@0: % matthiasm@0: % Description: matthiasm@0: % - applies an FIR filter matthiasm@0: % - takes care that f_audio and firfilter use same sampling frequency matthiasm@0: % - set samplingFreqFilter = samplingFreq if you do not which to resample matthiasm@0: % the filter matthiasm@0: % - predefines IR: 'GreatHall1','Classroom1','Octagon1', matthiasm@0: % 'GoogleNexusOneFrontSpeaker','GoogleNexusOneFrontMic' matthiasm@0: % 'VinylPlayer1960' matthiasm@0: % matthiasm@0: % Input: matthiasm@0: % f_audio - audio signal \in [-1,1]^{NxC} with C being the number of matthiasm@0: % channels matthiasm@0: % samplingFreq - sampling frequency of f_audio matthiasm@0: % timepositions_beforeDegr - some degradations delay the input signal. If matthiasm@0: % some points in time are given via this matthiasm@0: % parameter, timepositions_afterDegr will matthiasm@0: % return the corresponding positions in the matthiasm@0: % output. Set to [] if unavailable. Set f_audio matthiasm@0: % and samplingFreq to [] to compute only matthiasm@0: % timepositions_afterDegr. matthiasm@0: % matthiasm@0: % Input (optional): parameter matthiasm@0: % .presetIR = 'GreatHall1'; - load IR from a ADT preset matthiasm@0: % .impulseResponse = []; - IR to apply matthiasm@0: % .impulseResponseSampFreq = 0; - sampling rate of IR matthiasm@0: % .normalizeOutputAudio = 1 - normalize output matthiasm@0: % .normalizeImpulseResponse = 0 - l1-normalize firfilter before matthiasm@0: % application matthiasm@0: % matthiasm@0: % Output: matthiasm@0: % f_audio_out - audio signal \in [-1,1]^{NxC} with C being the number of matthiasm@0: % channels matthiasm@0: % timepositions_afterDegr - time positions corresponding to timepositions_beforeDegr matthiasm@0: % matthiasm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matthiasm@0: matthiasm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matthiasm@0: % Audio Degradation Toolbox matthiasm@0: % matthiasm@0: % Centre for Digital Music, Queen Mary University of London. matthiasm@0: % This file copyright 2013 Sebastian Ewert, Matthias Mauch and QMUL. matthiasm@0: % matthiasm@0: % This program is free software; you can redistribute it and/or matthiasm@0: % modify it under the terms of the GNU General Public License as matthiasm@0: % published by the Free Software Foundation; either version 2 of the matthiasm@0: % License, or (at your option) any later version. See the file matthiasm@0: % COPYING included with this distribution for more information. matthiasm@0: % matthiasm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matthiasm@0: matthiasm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matthiasm@0: % Check parameters matthiasm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matthiasm@0: if nargin<4 matthiasm@0: parameter=[]; matthiasm@0: end matthiasm@0: if nargin<3 matthiasm@0: timepositions_beforeDegr=[]; matthiasm@0: end matthiasm@0: if nargin<2 matthiasm@0: error('Please specify input data'); matthiasm@0: end matthiasm@0: if isfield(parameter,'loadInternalIR')==0 matthiasm@0: parameter.loadInternalIR = 1; matthiasm@0: end matthiasm@0: if isfield(parameter,'internalIR')==0 matthiasm@0: parameter.internalIR = 'GreatHall1'; matthiasm@0: end matthiasm@0: if isfield(parameter,'normalizeOutputAudio')==0 matthiasm@0: parameter.normalizeOutputAudio = 1; matthiasm@0: end matthiasm@0: if isfield(parameter,'normalizeImpulseResponse')==0 matthiasm@0: parameter.normalizeImpulseResponse = 0; matthiasm@0: end matthiasm@0: if isfield(parameter,'impulseResponse')==0 matthiasm@0: parameter.impulseResponse = []; matthiasm@0: end matthiasm@0: if isfield(parameter,'impulseResponseSampFreq')==0 matthiasm@0: parameter.impulseResponseSampFreq = 0; matthiasm@0: end matthiasm@0: matthiasm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matthiasm@0: % Main program matthiasm@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matthiasm@0: matthiasm@0: if parameter.loadInternalIR matthiasm@0: % load IR included in toolbox matthiasm@0: matthiasm@0: fullFilenameMfile = mfilename('fullpath'); matthiasm@0: [pathstr,name,ext] = fileparts(fullFilenameMfile); matthiasm@0: dirRootIRs = fullfile(pathstr,'degradationData'); matthiasm@0: matthiasm@0: names_internalIR = {'GreatHall1','Classroom1','Octagon1','GoogleNexusOneFrontSpeaker','GoogleNexusOneFrontMic','VinylPlayer1960'}; matthiasm@0: indexIR = find(strcmpi(names_internalIR,parameter.internalIR), 1); matthiasm@0: if isempty(indexIR) matthiasm@0: error('Please specify a valid preset') matthiasm@0: end matthiasm@0: matthiasm@0: switch indexIR matthiasm@0: case 1 matthiasm@0: file = fullfile(dirRootIRs,'RoomResponses/GreatHall_Omni_x06y06.wav'); matthiasm@0: case 2 matthiasm@0: file = fullfile(dirRootIRs,'RoomResponses/Classroom_Omni_30x20y.wav'); matthiasm@0: case 3 matthiasm@0: file = fullfile(dirRootIRs,'RoomResponses/Octagon_Omni_x06y06.wav'); matthiasm@0: case 4 matthiasm@0: file = fullfile(dirRootIRs,'PhoneResponses/IR_GoogleNexusOneFrontSpeaker.wav'); matthiasm@0: case 5 matthiasm@0: file = fullfile(dirRootIRs,'PhoneResponses/IR_GoogleNexusOneFrontMic.wav'); matthiasm@0: case 6 matthiasm@0: file = fullfile(dirRootIRs,'VinylSim/ImpulseReponseVinylPlayer1960_smoothed.wav'); matthiasm@0: end matthiasm@0: [parameter.impulseResponse,parameter.impulseResponseSampFreq] = wavread(file); matthiasm@0: end matthiasm@0: matthiasm@0: fs = samplingFreq; matthiasm@0: h = parameter.impulseResponse; matthiasm@0: fs_h = parameter.impulseResponseSampFreq; matthiasm@0: h_org = h; matthiasm@0: matthiasm@0: f_audio_out = []; matthiasm@0: if ~isempty(f_audio) matthiasm@0: matthiasm@0: if (size(h,2) > 1) && (size(h,1) ~= 1) matthiasm@0: error('Multichannel impulse responses are not supported'); matthiasm@0: end matthiasm@0: matthiasm@0: if fs ~= fs_h matthiasm@0: h = resample(h,fs,fs_h); matthiasm@0: end matthiasm@0: matthiasm@0: if parameter.normalizeImpulseResponse matthiasm@0: h = h / max(abs(h)); matthiasm@0: end matthiasm@0: matthiasm@0: f_audio_out = fftfilt(h,f_audio); matthiasm@0: matthiasm@0: if parameter.normalizeOutputAudio matthiasm@0: f_audio_out = adthelper_normalizeAudio(f_audio_out, samplingFreq); matthiasm@0: end matthiasm@0: end matthiasm@0: matthiasm@0: % This degradation does impose a delay matthiasm@0: timepositions_afterDegr = []; matthiasm@0: if ~isempty(timepositions_beforeDegr) matthiasm@0: % approximation via group delay matthiasm@0: matthiasm@0: [Gd,W] = mygrpdelay(h_org,1,fs_h,2048); matthiasm@0: %[Gd,W] = grpdelay(h_org,1,1024); % Matlab's own group delay function. Failt for some filters considerably. matthiasm@0: matthiasm@0: %figure; matthiasm@0: %plot(W,Gd/samplingFreqFilter) matthiasm@0: matthiasm@0: averageOfGroupDelays = mean(Gd); matthiasm@0: timeOffset_sec = averageOfGroupDelays / fs_h; matthiasm@0: matthiasm@0: timepositions_afterDegr = timepositions_beforeDegr + timeOffset_sec; matthiasm@0: end matthiasm@0: matthiasm@0: end matthiasm@0: matthiasm@0: function [gd,w] = mygrpdelay(b,a,Fs,nfft) matthiasm@0: % see also https://ccrma.stanford.edu/~jos/fp/Group_Delay_Computation_grpdelay_m.html matthiasm@0: matthiasm@0: b = b(:); matthiasm@0: a = a(:); matthiasm@0: matthiasm@0: w=Fs*[0:nfft-1]/nfft; matthiasm@0: oa = length(a)-1; % order of a(z) matthiasm@0: oc = oa + length(b)-1; % order of c(z) matthiasm@0: c = conv(b,fliplr(a)); % c(z) = b(z)*a(1/z)*z^(-oa) matthiasm@0: cr = c.*[0:oc]'; % derivative of c wrt 1/z matthiasm@0: num = fft(cr,nfft); matthiasm@0: den = fft(c,nfft); matthiasm@0: minmag = 10*eps; matthiasm@0: polebins = find(abs(den)