Update doa_music.m

Author: WenzheLiu-Speech

ssl_tools/doa_music.m

@@ -1,52 +1,34 @@
 function specGlobal = doa_music(x,Param,nsrc)
 if(size(x,2)<2)
-    error('ERROR[MUSIC]:信号通道数必须大于等于2');
+    error('ERROR[MUSIC]:信号通道数必须大于等于2');
 end
 %% STFT
 X = ssl_stft(x.',Param.window,Param.noverlap,Param.nfft,Param.fs);%nbin,nfram,nchan
 X = X(2:end,:,:);
 X = X(Param.freqBins,:,:);
 [nbin,~,nmic] = size(X);
 %% MUSIC
-% power = zeros(nbin, length(Param.azimuth)*length(Param.elevation));
-% for ibin = 1:nbin % 对于每个频点 
-%     Rxx = (transpose(squeeze(X(ibin,:,:)))*conj(squeeze(X(ibin,:,:))));% 自相关矩阵
-%     [U,~,~] = svd(Rxx);    % SVD分解   Rxx = U * S * U^H
-%     En = U(:,nsrc+1:end);  % 噪声子空间
-%     fprintf('%d\n',ibin)
-%     idx = 1;
-%     for iel = 1 : length(Param.elevation)%1 : length(Param.azimuth)
-%         for iaz = 1 : length(Param.azimuth)%1 : length(Param.elevation)
-%             % v为声源方向的单位向量，见鄢社锋书，需要注意的是，本程序的方位角与俯仰角与书中定义不同
-%             v = [cosd(Param.elevation(iel))*cosd(Param.azimuth(iaz));cosd(Param.elevation(iel))*sind(Param.azimuth(iaz));sind(Param.elevation(iel))];% 3 x 1
-%             tau = v'*(Param.micPos-repmat(Param.micPos(:,1),[1,nmic]))./Param.c; % 1 * nmic  参考麦克为1：Param.micPos(:,1)         
-%             a = exp(1i*2*pi*Param.f(ibin).*transpose(tau));%nmic x 1     SV = exp(-2*1i*pi*tau*Param.f.');  % nmic x nbin
-%             power(ibin,idx) = 1./(sum(abs( ctranspose(a) * En * ctranspose(En) * a )));
-%             idx = idx+1;
-%         end
-%     end
-% end
-% linspace包含端点，保证插值时不会出现NaN
+% linspace包含端点，保证插值时不会出现NaN
 aziGrid = linspace(Param.azimuth(1),Param.azimuth(end),round((Param.azimuth(end)-Param.azimuth(1))/Param.alphaRes)+1);
 eleGrid = linspace(Param.elevation(1),Param.elevation(end),round((Param.elevation(end)-Param.elevation(1))/Param.alphaRes)+1);
 power = zeros(nbin, length(aziGrid), length(eleGrid));
 
-for ibin = 1:nbin % 对于每个频点 
-    Rxx = (transpose(squeeze(X(ibin,:,:)))*conj(squeeze(X(ibin,:,:))));% 自相关矩阵
-    [U,~,~] = svd(Rxx);    % SVD分解   Rxx = U * S * U^H
-    En = U(:,nsrc+1:end);  % 噪声子空间
+for ibin = 1:nbin % 对于每个频点 
+    Rxx = (transpose(squeeze(X(ibin,:,:)))*conj(squeeze(X(ibin,:,:))));% 自相关矩阵
+    [U,~,~] = svd(Rxx);    % SVD分解   Rxx = U * S * U^H
+    En = U(:,nsrc+1:end);  % 噪声子空间
     fprintf('%d\n',ibin)
     for iaz = 1 :length(aziGrid)
         for iel = 1 :length(eleGrid)
             v = [cosd(eleGrid(iel))*cosd(aziGrid(iaz));cosd(eleGrid(iel))*sind(aziGrid(iaz));sind(eleGrid(iel))];% 3 x 1
-            tau = v'*(Param.micPos-repmat(Param.micPos(:,1),[1,nmic]))./Param.c; % 1 * nmic  参考麦克为1：Param.micPos(:,1)         
+            tau = v'*(Param.micPos-repmat(Param.micPos(:,1),[1,nmic]))./Param.c; % 1 * nmic  参考麦克为1：Param.micPos(:,1)         
             a = exp(1i*2*pi*Param.f(ibin).*transpose(tau));%nmic x 1     SV = exp(-2*1i*pi*tau*Param.f.');  % nmic x nbin
             power(ibin,iaz,iel) = 1./(sum(abs( ctranspose(a) * En * ctranspose(En) * a )));
         end
     end
 end
 
-% 对所有频率的空间谱加在一起:
+% 对所有频率的空间谱加在一起:
 spec = squeeze(sum(power,1)); %nAzi x nEle
 [az,el]=meshgrid(aziGrid,eleGrid);
 [azi,eli]=meshgrid(Param.azimuth,Param.elevation);
@@ -71,4 +53,4 @@
     X(:,:,ichan) = spectrogram(x(ichan,:),window,noverlap,nfft,fs); 
 end
 
-end
+end