Skip to content

GitLab

Commit cd83129b authored by Marco Berzborn's avatar Marco Berzborn
Browse files

fixed phase conventions for real valued spherical harmonic basis 

functions, added input parser to ita_sph_base
- Phase convention of the real valued SH basis now matches the AMBIX
convention (used in RAVEN)
- Deprecated SH basis functions removed
parent f30f2275
Hide whitespace changes
Inline Side-by-side
Showing with 183 additions and 336 deletions
applications/Binaural-HRTF/HRTF_class/@itaHRTF/interp.m
View file @ cd83129b
......@@ -106,7 +106,7 @@ ear_d = [-0.07 0.07];
W = sparse(diag(w)); % diagonal matrix containing weights
D = sparse(diag(dweights_rep)); % decomposition order-dependent Tikhonov regularization
Y = ita_sph_base(this.dirCoord,Nmax,'orthonormal',false); % calculate real-valued SHs using the measurement grid
Y = ita_sph_base(this.dirCoord,Nmax,'real'); % calculate real-valued SHs using the measurement grid
%% Calculate HRTF data for field points
if Nmax > 25
......@@ -131,7 +131,7 @@ for ear=1:2
% % calculate weighted SH coefficients using a decomposition order-dependent Tikhonov regularization
%
% freqData_temp = data.freqData;
% Y = ita_sph_base(data.channelCoordinates,Nmax,'orthonormal',false);
% Y = ita_sph_base(data.channelCoordinates,Nmax,'real');
a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.';
%a0 = (Y.'*W*Y + epsilon*D) \ Y.' *
......@@ -142,10 +142,10 @@ for ear=1:2
% calculate range-extrapolated HRTFs
a1 = a0 .* hankel_rep.';
Yest = ita_sph_base(field,N,'orthonormal',false); % use real-valued SH's
Yest = ita_sph_base(field,N,'real'); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a1).'; % interpolated + range-extrapolated HRTFs
else
Yest = ita_sph_base(field,Nmax,'orthonormal',false); % use real-valued SH's
Yest = ita_sph_base(field,Nmax,'real'); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a0).'; % interpolated HRTFs
end
end
......
applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m
View file @ cd83129b
......@@ -1110,7 +1110,7 @@ classdef itaHRTF < itaAudio
W = diag(vWeights); % diagonal matrix containing weights
D = diag(regweights_rep); % decomposition order-dependent Tikhonov regularization
Y = ita_sph_base(this.dirCoord,Nmeas,'orthonormal',false); % calculate real-valued SHs using the measurement grid (high SH-order)
Y = ita_sph_base(this.dirCoord,Nmeas,'real'); % calculate real-valued SHs using the measurement grid (high SH-order)
%% Calculate spatially smoothed HRTF data set
hrtf_smoo_wo_ITD = zeros(this.nBins,2*this.dirCoord.nPoints); % init.: columns: LRLRLR...
......@@ -1302,11 +1302,11 @@ classdef itaHRTF < itaAudio
earSidePlot = sArgs.earSide;
if numel(phiC_deg)>1,
xData = phiC_deg;
strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) ''];
strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) '�'];
strXlabel = '\phi in Degree';
else
xData = thetaC_deg;
strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) ''];
strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) '�'];
strXlabel = '\theta in Degree';
end
......@@ -1394,11 +1394,11 @@ classdef itaHRTF < itaAudio
earSidePlot = sArgs.earSide;
if numel(phiC_deg)>1,
xData = phiC_deg;
strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) ''];
strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) '�'];
strXlabel = '\phi in Degree';
else
xData = thetaC_deg;
strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) ''];
strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) '�'];
strXlabel = '\theta in Degree';
end
......
applications/SpatialAudio/SourceTranslation/m/calculateHandT_slideSource.m
View file @ cd83129b
......@@ -47,7 +47,7 @@ coordinate.cart = loc;
grid = params.display.grid + coordinate;
% Spherical Harmonics: num_angels x (Nc+1)^2.
Ynm = ita_sph_base(grid, N, 'Williams', true);
Ynm = ita_sph_base(grid, N);
% prepare hankel functions for translated source
for k_ind = 1 : K
......
applications/SpatialAudio/SourceTranslation/m/initParams.m
View file @ cd83129b
......@@ -168,8 +168,7 @@ params.display.grid.r = params.array.r;
% Spherical Harmonics: num_angels x (Nc+1)^2.
params.display.Ynm = ita_sph_base(params.display.grid,...
params.source.N,...
'Williams', true);
params.source.N);
% Build up a reverse grid.
th = params.display.grid.theta;
......@@ -207,8 +206,7 @@ params.array.mics = params.array.grid.nPoints; % Number of microphones
% Spherical Harmonics: num_angels x (Narray+1)^2.
params.array.Ynm = ita_sph_base(params.array.grid,...
params.array.N,...
'Williams', true);
params.array.N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters that are relative to the center
......
applications/SpatialAudio/SourceTranslation/m/simulateHandT_firstTranslation.m
View file @ cd83129b
......@@ -53,7 +53,7 @@ for k_ind = 1 : K
grid = params.display.grid - coordinate;
% Spherical Harmonics: num_angels x (Nc+1)^2.
Ynm = ita_sph_base(grid, N, 'Williams', true);
Ynm = ita_sph_base(grid, N);
hn_disPoints_x_N = ita_sph_besselh([0:N], 1,...
kr_disPoints_x_freqs(:,k_ind));
......
applications/SpatialAudio/ita_decodeAmbisonics.m
View file @ cd83129b
......@@ -57,7 +57,7 @@ for k=1:numel(weights)
Bformat(:,k)=weights(k).*Bformat(:,k); %Apply weighting factors
end
% SH and inversion
Y = ita_sph_base(LS,TheOrder,'williams',false); % generate basefunctions
Y = ita_sph_base(LS, TheOrder, 'real'); % generate basefunctions
Yinv=pinv(Y); % calculate Pseudoinverse
......
applications/SpatialAudio/ita_encodeAmbisonics.m
View file @ cd83129b
......@@ -17,7 +17,7 @@ opts.audioSource=''; %Signal of the source (itaAudio or filename)
opts = ita_parse_arguments(opts, varargin);
% Encode Source
Bformat = ita_sph_base(sourcePos, opts.order, 'Williams', false);
Bformat = ita_sph_base(sourcePos, opts.order, 'real');
if ~isempty(opts.audioSource)
if isa(opts.audioSource,'char')
......
applications/SphericalHarmonics/itaSamplingSphReal.m
View file @ cd83129b
......@@ -14,7 +14,7 @@ classdef itaSamplingSphReal < itaSamplingSph
end
methods(Access = protected)
function this = set_base(this, nmax)
this.Y = ita_sph_base(this,nmax,'Williams',false);
this.Y = ita_sph_base(this, nmax, 'real');
ita_verbose_info(['Real valued spherical harmonics up to order ' num2str(nmax) ' calculated.']);
end
end
......
applications/SphericalHarmonics/ita_sph_base.m
View file @ cd83129b
function Y = ita_sph_base(s, nmax, type, complex)
function Y = ita_sph_base(varargin)
%ITA_SPH_BASE - creates spherical harmonics (SH) base functions
% function Y = ita_sph_base(s, nmax, type)
% function Y = ita_sph_base(sampling, Nmax, options)
% Y is a matrix with dimensions [nr_points x nr_coefs]
% calculates matrix with spherical harmonic base functions
% calculates matrix with spherical harmonic basis functions
% for the grid given in theta and phi
% give type of normalization
%
%
% the definition was taken from:
......@@ -13,49 +12,51 @@ function Y = ita_sph_base(s, nmax, type, complex)
%
% This definition includes the Condon-Shotley phase term (-1)^m
%
% Martin Pollow (mpo@akustik.rwth-aachen.de)
% Institute of Technical Acoustics, RWTH Aachen, Germany
% 03.09.2008
% Syntax:
% Y = ita_sph_base(sampling, Nmax, options)
%
% Options (default):
% 'norm' ('orthonormal') : Normalization type
% 'real' (false) : Return real valued SH
%
% Example:
% Y = ita_sph_base(sampling, Nmax, 'norm', 'Williams')
%
% See also:
% ita_toolbox_gui, ita_read, ita_write, ita_generate
%
% Reference page in Help browser
% <a href="matlab:doc ita_sph_base">doc ita_sph_base</a>
%
% Deleted and new implementation
% (careful, history seems to be lost, but checkout of old version is possible):
% Bruno Masiero (bma@akustik.rwth-aachen.de)
% Institute of Technical Acoustics, RWTH Aachen, Germany
% 12.04.2011
% <ITA-Toolbox>
% This file is part of the application SphericalHarmonics for the ITA-Toolbox. All rights reserved.
% You can find the license for this m-file in the application folder.
% 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>
%
%
% Martin Pollow (mpo@akustik.rwth-aachen.de)
% Institute of Technical Acoustics, RWTH Aachen, Germany
% 03.09.2008
% check input
if nargin < 4
% Use complex version of the spherical harmonics
complex = true;
end
if nargin < 3
type = 'Williams';
% disp(['using default normalization: ' type]);
end
sArgs = struct('pos1_sampling', 'itaCoordinates', ...
'pos2_Nmax', 'int', ...
'norm', 'orthonormal', ...
'real', false);
[sampling, Nmax, sArgs] = ita_parse_arguments(sArgs, varargin);
%check for invalid sampling
if s.nPoints < 1
if sampling.nPoints < 1
Y = [];
ita_verbose_info('The sampling needs to consist of at least one point.', 0)
return
end
% vectorize grid angles
theta = s.theta;
phi = s.phi;
theta = sampling.theta;
phi = sampling.phi;
nr_points = numel(theta);
if nr_points ~= numel(phi)
error('theta and phi must have same number of points');
end
nr_coefs = (nmax+1)^2;
nm = 1:nr_coefs;
nm = 1:(Nmax+1)^2;
[~,m] = ita_sph_linear2degreeorder(nm);
exp_term = exp(1i*phi*m);
......@@ -63,16 +64,16 @@ exp_term = exp(1i*phi*m);
% calculate a matrix containing the associated legendre functions
% function Pnm = ass_legendre_func
% calculate the matrix of associated Legrendre functions
Pnm = zeros(nr_points, nr_coefs);
Pnm = zeros(sampling.nPoints, (Nmax+1)^2);
for ind = 0:nmax
for ind = 0:Nmax
% define the linear indices for the i'th degree
index_m_neg = ita_sph_degreeorder2linear(ind,-1:-1:-ind); % count in reverse order
index_m_pos = ita_sph_degreeorder2linear(ind,1:ind);
index_m_pos0 = ita_sph_degreeorder2linear(ind,0:ind);
% define positive orders with correct normalization
switch lower(type)
switch lower(sArgs.norm)
case {'orthonormal','williams'}
% the Pnm's used here are the Pnm's from Williams multiplied with the
% orthonormality factor sqrt((2n+1)./(4*pi).*(n-m)! ./ (n+m)!)
......@@ -85,7 +86,10 @@ for ind = 0:nmax
case {'schmidt','sch','semi-normalized'}
Pnm(:,index_m_pos0) = ...
bsxfun(@times,(-1).^(0:ind)*sqrt(2),legendre(ind,cos(theta.'),'sch').');
bsxfun(@times,(-1).^(0:ind),legendre(ind,cos(theta.'),'sch').');
case {'ambix'}
Pnm(:,index_m_pos0) = ...
bsxfun(@times,(-1).^(0:ind)/sqrt(8*pi),legendre(ind,cos(theta.'),'sch').');
otherwise
error('Wow! I do not know this normalization!')
......@@ -103,22 +107,12 @@ end
% compose the spherical harmonic base functions
Y = Pnm .* exp_term;
if ~complex
% now the conversion is done outside of this function
Y = ita_sph_complex2real(Y')';
% below the old code:
% for ind = 1:nmax
% % define the linear indices for the i'th degree
% index_m_neg = ita_sph_degreeorder2linear(ind,-1:-1:-ind); % count in reverse order
% index_m_pos = ita_sph_degreeorder2linear(ind,1:ind);
%
% for m = 1:length(index_m_neg)
% C = ((-1)^m*Y(:,index_m_pos(m)) + Y(:,index_m_neg(m)))/sqrt(2);
% S = ((-1)^m*Y(:,index_m_neg(m)) - Y(:,index_m_pos(m)))/sqrt(2)/1i;
%
% Y(:,ita_sph_degreeorder2linear(ind,m)) = C;
% Y(:,ita_sph_degreeorder2linear(ind,-m)) = S;
% end
% end
if strcmp(sArgs.norm, 'ambix')
% multiply by sqrt(2-delta_m)
mask = ~(m | zeros(size(m)));
Y(:,mask) = Y(:,mask) * sqrt(2);
end
if sArgs.real
Y = ita_sph_complex2real(Y.').';
end
\ No newline at end of file
applications/SphericalHarmonics/ita_sph_complex2real.m
View file @ cd83129b
function SH = ita_sph_complex2real(SH)
% converts a complex valued base or SH vector into its real base / SH vectors
function SH = ita_sph_complex2real(varargin)
%ITA_SPH_COMPLEX2REAL - Transforms from complex to valued basis functions
% This function transfroms spherical harmonic coefficients with complex
% valued basis functions to their corresponding real valued
% representation. The sign convention of the real valued basis functions
% can be chosen. The definition in 'Williams - Fourier Acoustics' is
% assumed for the complex basis functions (including the Condon-Shotley
% phase). The phase conventions for the real valued spherical harmonics
% can be 'raven' (which equals the ambix phase convention as defined in
% "Nachbar - AMBIX: A Suggested Ambisonics Format (revised by Zotter)")
% or 'zotter' as defined in "Zotter - Analysis and synthesis of
% sound-radiation with spherical arrays"
%
% Syntax:
% SH_real = ita_sph_complex2real(SH_cplx, options)
%
% Options (default):
% 'phase' ('ambix') : Phase definitions as in the AmbiX format
%
% Example:
% SH_real = ita_sph_complex2real(SH_cplx)
%
% See also:
% ita_toolbox_gui, ita_read, ita_write, ita_generate
%
% Reference page in Help browser
% <a href="matlab:doc ita_sph_complex2real">doc ita_sph_complex2real</a>
% <ITA-Toolbox>
% This file is part of the application SphericalHarmonics for the ITA-Toolbox. All rights reserved.
% You can find the license for this m-file in the application folder.
% 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>
% SH can be the vector/matrix of coefficients or the maximum order to
% calculate a matrix.
if isnumeric(SH)
if numel(SH) == 1
% return a matrix
SH = ita_sph_complex2real_mat(SH);
else
% assume this as SH coefs
SH = ita_sph_complex2real_coefs(SH);
end
else
error('please check syntax')
end
end
function SH = ita_sph_complex2real_coefs(SH)
% assume the SH dimension as first dimension
sizeY = size(SH);
% Rewrite, original author unknown
%
% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de
% Created: 12-Jun-2017
% check for max 2 dimension
if numel(sizeY) > 2
reshape(SH, sizeY(1), []);
end
sArgs = struct('pos1_SH', 'double', ...
'phase', 'ambix');
[SH, sArgs] = ita_parse_arguments(sArgs, varargin);
% check for SH in 1st dimension
if sizeY(1) == 1 && sizeY(2) > 1
SH = SH.';
isTransposed = true;
else
isTransposed = false;
if ~isnatural(sqrt(size(SH, 1)) - 1)
ita_verbose_info('The dimensions of the input data do not match.', 0)
SH = [];
return
end
nmax = sqrt(size(SH,1)) - 1;
if nmax ~= round(nmax)
error('ita_sph_complex2real something wrong with the number of coefficients (they do not fill up to a full order)');
end
Nmax = floor(sqrt(size(SH,1))-1);
for ind = 1:nmax
for ind = 1:Nmax
% define the linear indices for the i'th degree
index_m_neg = ita_sph_degreeorder2linear(ind,-1:-1:-ind); % count in reverse order
index_m_pos = ita_sph_degreeorder2linear(ind,1:ind);
for m = 1:length(index_m_neg)
cPos = SH(index_m_pos(m),:);
cNeg = SH(index_m_neg(m),:);
rPos = ((-1)^m * cPos + cNeg) / sqrt(2);
rNeg = ((-1)^m * cNeg - cPos) / sqrt(2) .* 1i;
switch sArgs.phase
case {'ambix', 'raven'}
rNeg = (cNeg - (-1)^m*cPos) / sqrt(2) .* 1i;
case 'zotter'
rNeg = (-cNeg + (-1)^m * cPos) / sqrt(2) .* 1i;
otherwise
ita_verbose_info('I do not know this phase convention. Aborting...', 0);
SH = [];
return
end
SH(ita_sph_degreeorder2linear(ind,m),:) = rPos;
SH(ita_sph_degreeorder2linear(ind,-m),:) = rNeg;
end
end
if isTransposed
SH = SH.';
end
% now bring to the old shape
SH = reshape(SH,sizeY);
end
function T = ita_sph_complex2real_mat(nmax)
if numel(nmax) > 1
error;
end
nSH = (nmax+1).^2;
lin = (1:nSH).';
multiplicator = 2*mod(lin,2) - 1;
[deg,order] = ita_sph_linear2degreeorder(lin);
linNeg = ita_sph_degreeorder2linear(deg,-order);
isPos = order > 0;
isNeg = order < 0;
T = sparse(eye(nSH));
T_orig = T;
T_swap = T .* diag(multiplicator);
T(lin(isPos),:) = (T_swap(lin(isPos),:) + T_orig(linNeg(isPos),:)) / sqrt(2);
T(lin(isNeg),:) = (T_swap(lin(isNeg),:) - T_orig(linNeg(isNeg),:)) / sqrt(2) .* 1i;
end
end
\ No newline at end of file
applications/SphericalHarmonics/ita_sph_mimo_error_simulation.m
View file @ cd83129b
......@@ -426,13 +426,12 @@ sArgs = struct('pos1_receiverCoords','itaCoordinates',...
'pos5_k','double',...
'r',[],...
'r_eq',1,...
'norm',false,...
'shType','complex');
'norm',false);
[receiverCoords,receiverNmax,sourceCoords,sourceNmax,k,sArgs] = ita_parse_arguments(sArgs,varargin);
[receiverLookDirection, sourceLookDirection,r] = array_orientation(receiverCoords,sourceCoords);
yReceiver = ita_sph_base(receiverLookDirection,receiverNmax,'Williams',strcmp(sArgs.shType,'complex'))';
ySource = ita_sph_base(sourceLookDirection,sourceNmax,'Williams',strcmp(sArgs.shType,'complex'))';
yReceiver = ita_sph_base(receiverLookDirection,receiverNmax)';
ySource = ita_sph_base(sourceLookDirection,sourceNmax)';
% avoid sArgs in parfor since it is a broadcast variable
if ~isempty(sArgs.r)
......
applications/SphericalHarmonics/ita_sph_real2complex.m
View file @ cd83129b
function SH = ita_sph_real2complex(SH)
% converts a real base or SH vector into its complex base / SH vectors
function SH = ita_sph_real2complex(varargin)
%ITA_SPH_REAL2COMPLEX - Transforms from real to complex valued basis functions
% This function transfroms spherical harmonic coefficients with real
% valued basis functions to their corresponding complex valued
% representation. The sign convention of the real valued basis functions
% can be chosen. The definition in 'Williams - Fourier Acoustics' is
% specified for the complex basis functions (including the Condon-Shotley
% phase). The phase conventions for the real valued spherical harmonics
% can be 'raven' (which equals the ambix phase convention as defined in
% "Nachbar - AMBIX: A Suggested Ambisonics Format (revised by Zotter)")
% or 'zotter' as defined in "Zotter - Analysis and synthesis of
% sound-radiation with spherical arrays"
%
% Syntax:
% SH_cplx = ita_sph_real2complex(SH_real, options)
%
% Options (default):
% 'phase' ('ambix') : Phase definitions as in the AmbiX format
%
% Example:
% SH_cplx = ita_sph_real2complex(SH_real)
%
% See also:
% ita_toolbox_gui, ita_read, ita_write, ita_generate
%
% Reference page in Help browser
% <a href="matlab:doc ita_sph_real2complex">doc ita_sph_real2complex</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>
% Rewrite, original author unknown
%
% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de
% Created: 12-Jun-2017
%
% <ITA-Toolbox>
......@@ -10,44 +46,20 @@ function SH = ita_sph_real2complex(SH)
% SH can be the vector/matrix of coefficients or the maximum order to
% calculate a matrix.
if isnumeric(SH)
if numel(SH) == 1
% return a matrix
SH = ita_sph_real2complex_mat(SH);
else
% assume this as SH coefs
SH = ita_sph_real2complex_coefs(SH);
end
else
error('please check syntax')
end
sArgs = struct('pos1_SH', 'double' ,...
'phase', 'ambix');
[SH, sArgs] = ita_parse_arguments(sArgs, varargin);
if ~isnatural(sqrt(size(SH, 1)) - 1)
ita_verbose_info('The dimensions of the input data do not match.', 0)
SH = [];
return
end
function SH = ita_sph_real2complex_coefs(SH)
Nmax = floor(sqrt(size(SH,1))-1);
% assume the SH dimension as first dimension
sizeY = size(SH);
% check for max 2 dimension
if numel(sizeY) > 2
reshape(SH, sizeY(1), []);
end
% check for SH in 1st dimension
if sizeY(1) == 1 && sizeY(2) > 1
SH = SH.';
isTransposed = true;
else
isTransposed = false;
end
nmax = sqrt(size(SH,1)) - 1;
if nmax ~= round(nmax)
error('ita_sph_complex2real something wrong with the number of coefficients (they do not fill up to a full order)');
end
for ind = 1:nmax
for ind = 1:Nmax
% define the linear indices for the i'th degree
index_m_neg = ita_sph_degreeorder2linear(ind,-1:-1:-ind); % count in reverse order
index_m_pos = ita_sph_degreeorder2linear(ind,1:ind);
......@@ -56,43 +68,19 @@ for ind = 1:nmax
rPos = SH(index_m_pos(m),:);
rNeg = SH(index_m_neg(m),:);
cPos = ((-1)^m .* rPos + 1i .* rNeg) ./ sqrt(2);
cNeg = (rPos - 1i .* (-1)^m .* rNeg) ./ sqrt(2);
switch lower(sArgs.phase)
case {'ambix','raven'}
cPos = (rPos + 1i.* rNeg) * (-1)^m / sqrt(2);
cNeg = (rPos - 1i * rNeg) ./ sqrt(2);
case {'zotter'}
cPos = (rPos - 1i.* rNeg) * (-1)^m / sqrt(2);
cNeg = (rPos + 1i * rNeg) ./ sqrt(2);
otherwise
ita_verbose_info('I do not know this phase convention. Aborting...', 0);
SH = [];
return
end
SH(ita_sph_degreeorder2linear(ind,m),:) = cPos;
SH(ita_sph_degreeorder2linear(ind,-m),:) = cNeg;
end
end
if isTransposed