fix for complex SHs

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

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

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

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;
 

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,

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;