% X_s = AKfractOctSmooth(x,[operation,fs,frac,mode,band,b,do_plot]);
%
% Applies fractional octave band smoothing on multichannel spectrum
% (for WELTI delivered as absolute or complex single-sided spectrum)
% (for SCHAERER delivered as impulse responses)
%
% To smooth a noisy spectrum try this
% x = AKnoise(2^12);
% X = AKboth2singleSidedSpectrum(fft(x));
% Y = AKfractOctSmooth(X,'amp',44100,3,true);
% % or
% Y = AKfractOctSmooth(x,'schaerer',44100,3,'amp','oct',[],true);
% y = ifft(AKsingle2bothSidedSpectrum(Y);
%
%
% INPUT operation - 'amp' for pure amplitude smoothing (zero phase),
% 'cmp' for seperate amplitude and phase smoothing,
% 'eqv' for amplitude smoothing plus copying phase of
% input signal,
% default = 'amp'
%
% output:
% X_s - smoothed single-sided spectrum with same length as the
% input spectrum
% input:
% x - single or multi channel SINGLE-SIDED SPECTRA
% fs - sampling frequency, default = 44100
% frac - fraction of octave for smoothing,
% default = 3 (third octave smoothing)
% do_plot - [0 .. 1 .. channel number] show comparative plot,
% [ fifth given as true/false or number of channel which should
% param ] be shown, default = false
%
% ----------------------
%
% INPUT operation - 'schaerer'
% - This algorithm implementation is quite expensive,
% especially for long input signals (> 4096 samples).
% Maybe apply signal truncation or windowing before!
%
% O U T P U T:
% X_s - smoothed single-sided spectrum with HALF EVEN length
% as the input spectrum
% input:
% x - single or multi channel SIMPULSE RESPONSES
% fs - sampling frequency, default = 44100
% frac - fraction of octave for smoothing,
% default = 3 (third octave smoothing)
% mode - 'amp' for amplitude,
% 'cpm' for complex,
% 'eqcmp' for equivalent complex (complex with amplitude
% compensation),
% default = 'amp'
% band - 'oct' for fractional octave,
% 'bark' for bark,
% 'ERB' for equivalent rectangular bandwidth,
% default = 'oct'
% b - coefficient for smoothing window
% W = b - (1-b)*cos(2*pi*k/N), for k=[0,length(x)-1],
% i.e. .5 (Hanning), .54 (Hamming), 1 (rectangular)],
% default = .5
% do_plot - [0 / 1 / channel number] show comparative plot,
% given as true/false or number(s) of channel which
% should be shown, default = false
%
% ----------------------------------------------------------------------- %
% Additional information to operational mode 'schaerer':
%
% Calculates the smoothed transfer function to an input [x]. Returns the
% smoothed transfer function X_s.
% The smoothed function is calculated by convolving the transfer function
% of H with a window function W. The width of the window increases with
% increasing frequency, [band] determines the frequency scale used
% (fraction of octave, bark or ERB). The convolution is either conducted
% with the magnitude transfer function or the complex transfer function. In
% the latter case, the window is also applied to the time domain,
% increasing in width with decreasing frequency. In the case of equivalent
% complex smoothing, X_s is a combination of the amplitude smoothed
% magnitude and the complex smoothed phase.
% [Hatziantoniou PD, Mourjopoulos JN (2000) Generalized Fractional-Octave
% Smoothing of Audio and Acoustic Responses]
% ----------------------------------------------------------------------- %
% v1 2009 Andreas Rotter, geeprombolo@gmail.com,
% Audio Communicatin Group, TU Berlin,
% using functions from Todd Welti and Andre Giese
% v2 01/2018 Hannes Helmholz, helmholz@campus.tu-berlin.de,
% Audio Communicatin Group, TU Berlin,
% rework and extention of Welti algorithm; improvement and
% smoothing fraction adjustment of Schaerer algorithm
% v3 02/2018 fixed behaviour for real spectrum inputs (e.g. AKp)
% v4 02/2018 calculation time warning and better string handling
% ----------------------------------------------------------------------- %

% AKtools
% Copyright (C) 2016 Audio Communication Group, Technical University Berlin
% Licensed under the EUPL, Version 1.1 or as soon they will be approved by
% the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at:
% http://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" basis,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

function X_s = AKfractOctSmooth(x,operation,fs,frac,mode,band,b,do_plot)

% an arbitrary value over which a Matlab warning gets shown in case the
% approximated calculation time for the Schaerer algorithm excals the given
% value
CALC_TIME_WARNING_SEC = 30;

if size(x,2) > size(x,1)
x = x.';
end

if nargin < 2 || isempty(operation)
operation = 'amp';
end
if nargin < 3 || isempty(fs)
fs = 44100;
end
if nargin < 4 || isempty(frac)
frac = 3;
end
if strcmpi(operation,'schaerer')
if nargin < 5 || isempty(mode)
mode = 'amp';
end
if nargin < 6 || isempty(band)
band = 'oct';
end
if nargin < 7 || isempty(b)
b = .5;
end
if nargin < 8
do_plot = 0;
end

switch lower(mode)
case 'amp'
fprintf('Schaerer -> amplitude smoothing');
case 'cmp'
fprintf('Schaerer -> complex smoothing');
case 'eqcmp'
fprintf('Schaerer -> equivalent complex smoothing');
otherwise
fprintf('\n');
error('AKfractOctSmooth:mode','Illegal mode.');
end
fprintf(' -> ');
switch lower(band)
case 'oct'
fprintf('bandwidth: 1/%g octave',frac);
case 'bark'
fprintf('bandwidth: critical band');
case 'erb'
fprintf('bandwidth: ERB');
otherwise
fprintf('\n');
error('AKfractOctSmooth:band','Illegal band.');
end
fprintf(' -> ...');
else
if nargin < 5
do_plot = 0;
else
do_plot = mode;
end
end

switch lower(operation)
% WELTI is a leftover from an old version
case {'amp','amplitude','abs','absolute','welti'}
operation = 'amp';
case {'cmp','complex','amp+pha','abs+pha','amp+ang','abs+ang'}
operation = 'cmp';
case {'eqv','equivalent','amp+cpy','abs+cpy'}
operation = 'eqv';
case 'schaerer'
% nothing to do here
otherwise
error('AKfractOctSmooth:operation','Illegal operation.');
end

% skip smoothing
if ~frac
if strcmpi(operation,'schaerer')
X_s = AKboth2singleSidedSpectrum(fft(x));
fprintf(' DONE.\n');
else
X_s = x;
end
return;
end

frac = 1/frac;

% ----------------------------------------------------------------------- %

if strcmpi(operation,'schaerer')
time = SchaererApproxCalcTime(size(x,1),size(x,2));
if time > CALC_TIME_WARNING_SEC
fprintf('\n');
warning('AKfractOctSmooth:schaerer',['The estimated '...
'calculation time for %d impulse response(s) of %d samples '...
'is about %d seconds.'],size(x,2),size(x,1),ceil(time));
fprintf(' -> ...');
end
X_s = AKfractOctSmoothSchaerer(fft(x),fs,frac,mode,band,b);

% ------------------------------------------------------------------- %
else
% frequency support vector
f = linspace(0,fs/2,size(x,1));

if isreal(x) && strcmpi(operation,'amp')
% direct smoothing, used e.g. by AKp
X_s = AKfractOctSmoothWelti(x,frac,f(end),f);

else
% seperate magnitude and phase
X_s_abs = AKfractOctSmoothWelti(abs(x).^2,frac,f(end),f);
switch lower(operation)
case 'amp'
X_s_angle = zeros(size(x));
case 'cmp'
X_s_angle = AKfractOctSmoothWelti(angle(x),frac,f(end),f);
% % One could do the smoothing on the unwrapped phase
% % which leads to a different but not more meaningful IR.
% % Panzer2004 ï¿½The Use of Continuous Phase for ..."
% % suggests this should work somehow though.
% % (H. Helmholz, 01/2018)
% X_s_angle = AKfractOctSmoothWelti(unwrap(angle(x)),frac,f(end),f);
case 'eqv'
X_s_angle = angle(x);
end
X_s = sqrt(X_s_abs) .* exp(1i*X_s_angle);
end
end

% ----------------------------------------------------------------------- %

% plotting
if do_plot
channels = do_plot; % also the channel index which should be plotted
if strcmpi(operation,'schaerer')
operation_plot = [mode,' ',band];
x_plot = x(:,channels);
else
operation_plot = operation;
x_plot = ifft(AKsingle2bothSidedSpectrum(x(:,channels),true));
end
x_s_plot = ifft(AKsingle2bothSidedSpectrum(X_s(:,channels),true));

AKf;
set(gcf,'Name',sprintf('%s smoothing (ch %s)',...
operation_plot,AKlistIntegers(channels)),'NumberTitle','off');
AKpMulti(x_plot,{'t2d','et2d';'m2d','p2d'},'c','k');
AKpMulti(x_s_plot,{'t2d','et2d';'m2d','p2d'},'c','r','lw',2);
AKtightenFigure;
end

if strcmpi(operation,'schaerer')
fprintf(' DONE.\n');
end

end

% ----------------------------------------------------------------------- %
function H_smooth = AKfractOctSmoothWelti(H,oct,f_stop_i,f)

if nargin >= 4
f_stop_i = round(f_stop_i/f(2))+1;
end

% cast H to double to avoid error (F. Brinkmann 7/2012)
H = cast(H,'double');

% fixed at 1 for now, doesn't work if not starting at 1
f_start_i = 1;
% number of sam ples for log warping, seems to work well
f_N = f_stop_i;
% this is a multiplicitive factor, not an even spacing
oct_spacing = 10^(log10(f_stop_i-f_start_i)/f_N);
% number of bins per xth octave
oct_N = oct*log10(2)/log10(oct_spacing);

% logarithmic warping
f_log = logspace(log10(f_start_i),log10(f_stop_i),f_N);
H_log = interp1(1:length(H),H,f_log,'spline');
H_log = reshape(H_log,size(H));
% don't change to SPLINE(), because of interpolation per column (not row)

% make it even length for convenience
oct_N_even = round(oct_N/2)*2;
win = gausswin(oct_N_even*2);

% minimze smoothing error at boundaries by extension
H_ext = ones(length(win),size(H_log,2),size(H_log,3));
H_log = [H_ext .* H_log(1);H_log;H_ext .* H_log(end)];

% convolve with smothing window
H_smooth = zeros(size(H_log));
for ch = 1:size(H_log,ndims(H_log))
H_smooth(:,ch) = fftfilt(win,H_log(:,ch),2*length(win))./sum(win);
end

% get rid of extra bins from adding lead and lag
% H_smooth = H_smooth(length(H_ext)+.5*length(win) + [1:length(f_log)],:,:); %#ok<NBRAK>
H_smooth = H_smooth(size(H_ext,1)+.5*size(win,1) + [1:length(f_log)],:,:); %EDITED ISAAC

% dewarping to linear
H_smooth = interp1(f_log,H_smooth,1:f_stop_i,'spline');
H_smooth = reshape(H_smooth,size(H));
% don't change to SPLINE(), because of interpolation per column (not row)

% transpose if needed to get same output as input (column/row)
if size(H,1)~=f_stop_i && size(H,2)~=1
H_smooth = H_smooth';
end

end

% ----------------------------------------------------------------------- %
function H_smooth = AKfractOctSmoothSchaerer(H,fs,frac,mode,band,b)

% reshape matrix to also handle [filterlen,1,channels] inputs
H_orig_size = size(H);
H = squeeze(H);

[filterlen,channels] = size(H);
filterlen_2 = floor(filterlen/2);

if any(strcmpi(band,{'bark','erb'}))
frac = 1/frac;
b = frac;
end

% generate vector with bandwidths according to the selected scale
d_f = fs/filterlen;
switch lower(band)
case 'oct'
% % original Schaerer implementation
% f_up = (10^(3/10)) ^ (.5*frac);
% f_low = (10^(-3/10)) ^ (.5*frac);
% bandwidth = (1:filterlen_2) * d_f*(f_up-f_low);

% changes see below (H. Helmholz 01/2018)
bandwidth = (1:filterlen_2) * d_f;
case 'bark'
bandwidth = 25+75*(1+1.4*((1:filterlen_2) * d_f/1000).^2) .^ .69;
case 'erb'
bandwidth = 24.7*(4.37*((1:filterlen_2) * d_f/1000) + 1);
end
% the sqrt(3*frac) correction factor was imcluded, to match the Welti
% algorithm smoothing behaviour (with b=.5 and band='oct')
% (H. Helmholz 01/2018)
m = ceil(.5 * bandwidth * sqrt(3*frac)/d_f);

% generate windowfunction
W = zeros(filterlen_2,filterlen);
for s = 1:filterlen_2
M = m(s);
W(s,1:M+1) = (b - (b-1)*cos((pi/M)*(0:M))) ./ (2*b*(M+1)-1);
W(s,filterlen-M+1:filterlen) = ...
(b - (b-1)*cos((pi/M)*((filterlen-M:filterlen-1) - filterlen))) ...
./ (2*b*(M+1)-1);
end

% seperate amplitde and phase
if any(strcmpi(mode,{'amp','eqcmp'}))
H_2 = abs(H).^2;

% generate circular convolution matrix
H_circ_2 = zeros(filterlen,filterlen,channels);
for s = 1:filterlen
H_circ_2(s,:,:) = circshift(H_2,s-1,1);
end

% circular convolution
H_smooth_circ_2 = zeros(filterlen_2,filterlen,channels);
for ch = 1:channels
H_smooth_circ_2(:,:,ch) = W * H_circ_2(:,:,ch);
end

% derive smoothed function
H_smooth_2 = zeros(filterlen,channels);
for s = 1:filterlen_2
H_smooth_2(s,:) = squeeze(H_smooth_circ_2(m(s),s,:));
if s > 1 % mirror spectrum
H_smooth_2(filterlen-s+2,:) = H_smooth_2(s,:);
end
end
end
if any(strcmpi(mode,{'cmp','eqcmp'}))
H_r = real(H);
H_im = imag(H);

% generate circular convolution matrix
H_circ_r = zeros(filterlen,filterlen,channels);
H_circ_im = H_circ_r;
for s = 1:filterlen
H_circ_r(s,:,:) = circshift(H_r,s-1,1);
H_circ_im(s,:,:) = circshift(H_im,s-1,1);
end

% circular convolution
H_smooth_circ_r = zeros(filterlen_2,filterlen,channels);
H_smooth_circ_im = H_smooth_circ_r;
for ch = 1:channels
H_smooth_circ_r(:,:,ch) = W * H_circ_r(:,:,ch);
H_smooth_circ_im(:,:,ch) = W * H_circ_im(:,:,ch);
end

% derive smoothed function
H_smooth_r = zeros(filterlen,channels);
H_smooth_im = H_smooth_r;
for s = 1:filterlen_2
H_smooth_r(s,:) = squeeze(H_smooth_circ_r(m(s),s,:));
H_smooth_im(s,:) = squeeze(H_smooth_circ_im(m(s),s,:));
if s > 1 % mirror spectrum
H_smooth_r(filterlen-s+2,:) = H_smooth_r(s,:);
H_smooth_im(filterlen-s+2,:) = H_smooth_im(s,:);
end
end
H_smooth_c = H_smooth_r + 1i*H_smooth_im;
end

% combine amplitude and phase
switch lower(mode)
case 'amp'
H_smooth = sqrt(H_smooth_2) .* exp(1i*angle(H));
case 'cmp'
H_smooth = H_smooth_c;
case 'eqcmp'
H_smooth = sqrt(H_smooth_2) .* exp(1i*angle(H_smooth_c));
end

% restore input dimensions
H_smooth = reshape(H_smooth,H_orig_size);
H_smooth = AKboth2singleSidedSpectrum(H_smooth);

end

function time_sec = SchaererApproxCalcTime(ir_len,n_irs)
% time approximation based on these values on a 3.5 GHz 6-Core machine:
%
% len = [1024,2048,4096,8192,16384];
% time = [0.518675,0.92,3.2,16.205843,129.771364];
% weights = [.1,.25,.5,1,1];
%
% which leads to the following coefficient by curve fitting:

a = 1.31e-11;
b = 3.083;
c = 1.052;

time_sec = n_irs * a * ir_len ^ b + c;

% This approximation will be not quite right for another machine of course.
% Also the scaling with an increasing number of IRs is not linear (despite
% being calculated like that) since some parts of the algorithm will be
% parallelized by Matlab already.

end