From 6eef482dfe8457f5dc13a2656595d9bd521ea0b8 Mon Sep 17 00:00:00 2001 From: Thomas Deppisch Date: Fri, 20 May 2022 10:57:28 +0200 Subject: [PATCH] fix for complex SHs --- dependencies/getSMAIRMatrix.m | 6 +++--- getEMagLs2Filters.m | 20 +++++++------------- getEMagLsFilters.m | 23 +++++++++++------------ getLsFilters.m | 6 +++--- getMagLsFilters.m | 32 ++++++++++++++++---------------- testEMagLs.m | 2 +- 6 files changed, 41 insertions(+), 48 deletions(-) diff --git a/dependencies/getSMAIRMatrix.m b/dependencies/getSMAIRMatrix.m index d1ac4cd..49a96c4 100644 --- a/dependencies/getSMAIRMatrix.m +++ b/dependencies/getSMAIRMatrix.m @@ -106,17 +106,17 @@ for ii = 1:numFreqs Bn = diag(sh_repToOrder(bnAll(ii,:).').'); - pMics(:,:,ii) = YMics' * Bn; + pMics(:,:,ii) = YMics.' * Bn; if ii > 1 && ii < numFreqs - pMics(:,:,end-ii+2) = YMics' * conj(Bn); + pMics(:,:,end-ii+2) = YMics.' * conj(Bn); end end pMics(isnan(pMics)) = 0; % remove singularity for f = 0 pN = zeros(numShsOut, numShsSimulation, nfft); - pinvYMicsLow = pinv(YMicsLow'); + pinvYMicsLow = pinv(YMicsLow.'); for ii = 1:nfft pN(:,:,ii) = pinvYMicsLow * pMics(:,:,ii); end diff --git a/getEMagLs2Filters.m b/getEMagLs2Filters.m index bffcc0b..17c31d0 100644 --- a/getEMagLs2Filters.m +++ b/getEMagLs2Filters.m @@ -37,7 +37,7 @@ numDirections = size(hL,2); nfft = max(2048,2*len); simulationOrder = 32; -YHi = getSH(simulationOrder, [hrirGridAziRad,hrirGridZenRad], shDefinition).'; +YHi = conj(getSH(simulationOrder, [hrirGridAziRad,hrirGridZenRad], shDefinition)); f_cut = 2000; % transition frequency f = linspace(0,fs/2,nfft/2+1); @@ -77,26 +77,26 @@ W_MLS_l = zeros(nfft, numMics); W_MLS_r = zeros(nfft, numMics); for k = 1:numPosFreqs % leave out first bin, here bn = 0 - pwGrid = smairMat(:,:,k) * conj(YHi); + pwGrid = smairMat(:,:,k) * YHi.'; [U,S,V] = svd(pwGrid.','econ'); s = diag(S); s = 1 ./ max(s, svdRegulConst * max(s)); % regularize regInvY = (V .* s.') * U'; - W_LS_l(k,:) = (regInvY * HL(k,:).').'; - W_LS_r(k,:) = (regInvY * HR(k,:).').'; + W_LS_l(k,:) = HL(k,:) * regInvY.'; + W_LS_r(k,:) = HR(k,:) * regInvY.'; % calculate negative frequencies (important in case of complex-valued % SHs) if k > 1 && k < numPosFreqs - pwGridNeg = smairMat(:,:,end-k+2) * conj(YHi); + pwGridNeg = smairMat(:,:,end-k+2) * YHi.'; [UNeg,SNeg,VNeg] = svd(pwGridNeg.','econ'); sNeg = diag(SNeg); sNeg = 1 ./ max(sNeg, svdRegulConst * max(sNeg)); % regularize regInvYNeg = (VNeg .* sNeg.') * UNeg'; - W_LS_l(end-k+2,:) = (regInvYNeg * conj(HL(k,:)).').'; - W_LS_r(end-k+2,:) = (regInvYNeg * conj(HR(k,:)).').'; + W_LS_l(end-k+2,:) = conj(HL(k,:)) * regInvYNeg.'; + W_LS_r(end-k+2,:) = conj(HR(k,:)) * regInvYNeg.'; end if k >= k_cut @@ -170,13 +170,7 @@ W_MLS_r = conj(HCorr(:,:,2)); end -W_MLS_l(1,:) = real(W_MLS_l(2,:)); % DC extension -W_MLS_r(1,:) = real(W_MLS_r(2,:)); -W_MLS_l(end,:) = real(W_MLS_l(end,:)); -W_MLS_r(end,:) = real(W_MLS_r(end,:)); -W_MLS_l = [W_MLS_l; flipud(conj(W_MLS_l(2:end-1,:)))]; wMlsL = ifft(W_MLS_l,nfft); -W_MLS_r = [W_MLS_r; flipud(conj(W_MLS_r(2:end-1,:)))]; wMlsR = ifft(W_MLS_r,nfft); % shorten, shift diff --git a/getEMagLsFilters.m b/getEMagLsFilters.m index ce2f5bc..f34086a 100644 --- a/getEMagLsFilters.m +++ b/getEMagLsFilters.m @@ -38,21 +38,20 @@ numDirections = size(hL,2); nfft = max(2048,2*len); simulationOrder = 32; -YHi = getSH(simulationOrder, [hrirGridAziRad,hrirGridZenRad], shDefinition).'; -YLo = YHi(1:numHarmonics,:); -pinvYLo = pinv(YLo'); +YHi = conj(getSH(simulationOrder, [hrirGridAziRad,hrirGridZenRad], shDefinition)); +YLo = YHi(:,1:numHarmonics); +pinvYLo = pinv(YLo.'); f_cut = 2000; % transition frequency f = linspace(0,fs/2,nfft/2+1); k_cut = round(f_cut/f(2) + 1); % zero pad and remove group delay (alternative to applying global phase delay later) -grpDL = grpdelay((pinvYLo(1,:) * hL.').', 1, f, fs); -grpDR = grpdelay((pinvYLo(1,:) * hR.').', 1, f, fs); +grpDL = grpdelay(hL * pinvYLo(:,1), 1, f, fs); +grpDR = grpdelay(hR * pinvYLo(:,1), 1, f, fs); hL = circshift([hL; zeros(nfft - size(hL, 1), size(hL, 2))], -round(median(grpDL))); hR = circshift([hR; zeros(nfft - size(hR, 1), size(hR, 2))], -round(median(grpDR))); - HL = fft(hL,nfft); HR = fft(hR,nfft); @@ -80,26 +79,26 @@ W_MLS_l = zeros(nfft, numHarmonics); W_MLS_r = zeros(nfft, numHarmonics); for k = 1:numPosFreqs % leave out first bin, here bn = 0 - pwGrid = smairMat(:,:,k) * conj(YHi); + pwGrid = smairMat(:,:,k) * YHi.'; [U,S,V] = svd(pwGrid.','econ'); s = diag(S); s = 1 ./ max(s, svdRegulConst * max(s)); % regularize regInvY = (V .* s.') * U'; - W_LS_l(k,:) = (regInvY * HL(k,:).').'; - W_LS_r(k,:) = (regInvY * HR(k,:).').'; + W_LS_l(k,:) = HL(k,:) * regInvY.'; + W_LS_r(k,:) = HR(k,:) * regInvY.'; % calculate negative frequencies (important in case of complex-valued % SHs) if k > 1 && k < numPosFreqs - pwGridNeg = smairMat(:,:,end-k+2) * conj(YHi); + pwGridNeg = smairMat(:,:,end-k+2) * YHi.'; [UNeg,SNeg,VNeg] = svd(pwGridNeg.','econ'); sNeg = diag(SNeg); sNeg = 1 ./ max(sNeg, svdRegulConst * max(sNeg)); % regularize regInvYNeg = (VNeg .* sNeg.') * UNeg'; - W_LS_l(end-k+2,:) = (regInvYNeg * conj(HL(k,:)).').'; - W_LS_r(end-k+2,:) = (regInvYNeg * conj(HR(k,:)).').'; + W_LS_l(end-k+2,:) = conj(HL(k,:)) * regInvYNeg.'; + W_LS_r(end-k+2,:) = conj(HR(k,:)) * regInvYNeg.'; end if k >= k_cut diff --git a/getLsFilters.m b/getLsFilters.m index 36d0ae1..e90b91c 100644 --- a/getLsFilters.m +++ b/getLsFilters.m @@ -22,8 +22,8 @@ shDefinition = 'complex'; % real or complex Y = getSH(order, [hrirGridAziRad,hrirGridZenRad], shDefinition); -pinvY = pinv(Y); +pinvY = pinv(Y'); -wLsL = hL * pinvY'; -wLsR = hR * pinvY'; +wLsL = hL * pinvY; +wLsR = hR * pinvY; diff --git a/getMagLsFilters.m b/getMagLsFilters.m index cd62228..713fd07 100644 --- a/getMagLsFilters.m +++ b/getMagLsFilters.m @@ -31,24 +31,24 @@ end nfft = max(2*len,2048); -Y = getSH(order, [hrirGridAziRad,hrirGridZenRad], shDefinition); -pinvY = pinv(Y); +Y = conj(getSH(order, [hrirGridAziRad,hrirGridZenRad], shDefinition)); +pinvY = pinv(Y.'); f_cut = 2000; % transition frequency f = linspace(0,fs/2,nfft/2+1); k_cut = round(f_cut/f(2) + 1); % zero pad and remove delay (alternative to applying global phase delay later) -grpDL = grpdelay(hL * pinvY(1,:)', 1, f, fs); -grpDR = grpdelay(hR * pinvY(1,:)', 1, f, fs); +grpDL = grpdelay(hL * pinvY(:,1), 1, f, fs); +grpDR = grpdelay(hR * pinvY(:,1), 1, f, fs); hL = circshift([hL; zeros(nfft - size(hL, 1), size(hL, 2))], -round(median(grpDL))); hR = circshift([hR; zeros(nfft - size(hR, 1), size(hR, 2))], -round(median(grpDR))); HL = fft(hL,nfft); HR = fft(hR,nfft); -w_LS_l = hL * pinvY'; -w_LS_r = hR * pinvY'; +w_LS_l = hL * pinvY; +w_LS_r = hR * pinvY; W_LS_l = fft(w_LS_l,nfft); W_LS_r = fft(w_LS_r,nfft); @@ -58,25 +58,25 @@ W_MLS_l = W_LS_l; W_MLS_r = W_LS_r; for k = k_cut:numPosFreqs-1 - phi_l = angle(W_MLS_l(k-1,:) * Y'); - W_MLS_l(k,:) = (abs(HL(k,:)) .* exp(1i * phi_l)) * pinvY'; + phi_l = angle(W_MLS_l(k-1,:) * Y.'); + W_MLS_l(k,:) = (abs(HL(k,:)) .* exp(1i * phi_l)) * pinvY; - phi_r = angle(W_MLS_r(k-1,:) * Y'); - W_MLS_r(k,:) = (abs(HR(k,:)) .* exp(1i * phi_r)) * pinvY'; + phi_r = angle(W_MLS_r(k-1,:) * Y.'); + W_MLS_r(k,:) = (abs(HR(k,:)) .* exp(1i * phi_r)) * pinvY; % calculate negative frequencies (important in case of complex-valued % SHs) - phi_l_n = angle(conj(W_MLS_l(k-1,:)) * Y'); - W_MLS_l(end-k+2,:) = (abs(HL(k,:)) .* exp(1i * phi_l_n)) * pinvY'; + phi_l_n = angle(conj(W_MLS_l(k-1,:)) * Y.'); + W_MLS_l(end-k+2,:) = (abs(HL(k,:)) .* exp(1i * phi_l_n)) * pinvY; - phi_r_n = angle(conj(W_MLS_r(k-1,:)) * Y'); - W_MLS_r(end-k+2,:) = (abs(HR(k,:)) .* exp(1i * phi_r_n)) * pinvY'; + phi_r_n = angle(conj(W_MLS_r(k-1,:)) * Y.'); + W_MLS_r(end-k+2,:) = (abs(HR(k,:)) .* exp(1i * phi_r_n)) * pinvY; end % Nyuist bin k = numPosFreqs; -phi_l = angle(W_MLS_l(k-1,:) * Y'); -W_MLS_l(k,:) = real(abs(HL(k,:)) .* exp(1i * phi_l)) * pinvY'; +phi_l = angle(W_MLS_l(k-1,:) * Y.'); +W_MLS_l(k,:) = real(abs(HL(k,:)) .* exp(1i * phi_l)) * pinvY; if applyDiffusenessConst % diffuseness constraint after Zaunschirm, Schoerkhuber, Hoeldrich, diff --git a/testEMagLs.m b/testEMagLs.m index c6bf318..890371d 100644 --- a/testEMagLs.m +++ b/testEMagLs.m @@ -42,7 +42,7 @@ %% sh transform and radial filter (for LS and conventional MagLS) E = getSH(shOrder, [micGridAziRad, micGridZenRad], 'complex'); -shRecording = smaRecording * pinv(conj(E)).'; +shRecording = smaRecording * pinv(E).'; % parameters for the radial filter params.order = shOrder;