Macierze mikrofonowe i głośnikowe

Autorzy: Daniel Król

Pliki: rozdzial_17.zip

Przykład 17.1

W przykładzie 17.1 (plik przyklad_17_1.m) przedstawiono program w języku MATLAB, który umożliwia wygenerowanie charakterystyk kierunkowych dla różnej liczby mikrofonów w macierzy na podstawie równania 17.1 oraz 17.4 (kod zakomentowany).
% Przykład 17.1
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

clear all; clc;                         % przygotowanie obliczeń
%------------------------ Parametry macierzy ------------------------------
v=340;                                  % prędkość dźwięku [m/s]
l=0.05;                                 % odległość mikrofonów [m]
N=8;                                    % liczba mikrofonów
f=1000;                                 % częstotliwość [Hz]
%--------------------------------------------------------------------------
m=1;
for theta=-pi:0.01:pi                   % zmiana kąta -pi do +pi
    sum=0;
    %re=0; im=0;
    for k=0:N-1                         % obliczenie czułości macierzy
        sum=sum+exp(i*2*pi*f*k*l*sin(theta)/v); % na podstawie wzoru (17.1)
        %re=re+cos(2*pi*f*k*l*sin(theta)/v);    % na podstawie wzoru (17.4)
        %im=im+sin(2*pi*f*k*l*sin(theta)/v);    % na podstawie wzoru (17.4)
    end
    y(m)=(abs(sum/N));                  % moduł/N - na podstawie wzoru (17.1)
    %y(m)=sqrt(re*re+im*im)/N;          % moduł/N - na podstawie wzoru (17.4)
    m=m+1;
end
%--------------------------------------------------------------------------
ph=linspace(-pi,pi,m-1);                % przestrzeń fazy dla wykresu
polar(ph,y,'b');                        % wykres biegunowy
xlabel(['N=',num2str(N),', l=',num2str(l*100),'cm', ', f=',num2str(f/1000),'kHz']);

Przykład 17.2

W przykładzie 17.2 (plik przyklad_17_2.m) przedstawiono program w języku MATLAB, który umożliwia wygenerowanie szerokopasmowej charakterystyki kierunkowej dla różnej liczby mikrofonóww macierzy na podstawie wzoru 17.1.
% Przykład 17.2
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

clear all; clc;                         % przygotowanie obliczeń
%------------------------ Parametry macierzy ------------------------------
v=340;                                  % prędkość dźwięku [m/s]
l=0.05;                                 % odległość mikrofonów [m]
N=8;                                    % liczba mikrofonów
fg=6000;                                % górna częstotliwość wykresu [Hz]
M=400;                                  % rozdzielczość wykresu MxM
%------------------------ Obliczenia czułości -----------------------------
for m=1:M                               % zmiana częstotliwości
    f(m)=m*fg/M;                        % wektor częstotliwości
    for p=1:M                           % zmiana fazy
        theta(p)=(pi*(p/M))-pi/2;       % wektor fazy
        sum=0;
        for k=0:N-1                     % obliczenie czułości macierzy
            sum=sum+exp(i*2*pi*f(m)*k*l*sin(theta(p))/v);
        end
        y(m,p)=abs(sum/N);              % moduł/N
    end
end
%----------------------------- Wykres -------------------------------------
fr=linspace(0,fg,M);                    % przestrzeń częstotliwości dla wykresu
ph=linspace(-90,90,M);                  % przestrzeń fazy dla wykresu
%surf(ph,fr,y);                         % wykres 3D
pcolor(ph,fr,y); colorbar;              % wykres 2D
shading interp;                         % cieniowanie wykresu
ylabel('Częstotliwość (Hz)'); xlabel('Faza (stopnie)'); zlabel('Amplituda');
title(['N=',num2str(N),', l=',num2str(l*100),'cm']); grid on;

Przykład 17.3

W przykładzie 17.3 (plik delay_sum.m) przedstawiono implementację w języku MATLAB algorytmu Delay-Sum.
% Przykład 17.3
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

function y=delay_sum(X, fp, l, theta)
%% y=delay_sum(X, fp, l, theta)
%
% X - wielokanałowy sygnał wejsciowy, fp - częstotliwość próbkowania
% l - odległość mikrofonów, theta - kąt wiązki
% y - sygnał wyjściowy
v=340;                                  % prędkość dźwięku [m/s]
alfa=(theta*pi)/180;                    % stopnie/radiany
d=(l*sin(alfa))/v;                      % opóźnienie [s]
N=size(X,2);                            % liczba kanałów
%--------------------------------------------------------------------------
for k=1:N
    shf(N-k+1)=floor((k-1)*abs(d)*fp);  % wektor opóźnień
end
if(d<0)                                 % transformacja dla
    shf=fliplr(shf);                    % ujemnego kąta wiązki
end
%--------------------------------------------------------------------------
M=size(X,1);                            % liczba próbek w kanale
y1=zeros(M,1);
for k=1:N
    y1(1+shf(k):M+shf(k),k)=X(1:M,k);   % opóźnienia
end
s=ones(1,N)/N;                          % przygotowanie wektora sumującego
y=y1*s';                                % suma opóźnionych kanałów / N

Przykład 17.4

W przykładzie 17.4 (plik filter_sum.m) przedstawiono implementację w języku MATLAB algorytmu Filter-Sum.
% Przykład 17.4
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

function y=filter_sum(X, fp, h, l, theta)
%% y=filter_sum(X, fp, h, l, theta)
%
% X - wielokanałowy sygnał wejściowy, fp - częstotliwość próbkowania
% h - wagi podfiltru, l - odległość mikrofonów, theta - kąt wiązki
% y - sygnał wyjściowy
v=340;                                  % prędkość dźwięku
alfa=(theta*pi)/180;                    % stopnie/radiany
d=(l*sin(alfa))/v;                      % opóźnienie [s]
N=size(X,2);                            % liczba kanałów
%------------------------------------------------------------
for k=1:N
    shf(N-k+1)=floor((k-1)*abs(d)*fp);  % wektor opóźnień
end
if(d<0)                                 % transformacja dla
    shf=fliplr(shf);                    % ujemnego kąta wiązki
end
%------------------------------------------------------------
ws=size(h,2);                           % długość podfiltru
dmax=floor(fp*(((k-1)*l)/v));           % max opóźnienie w macierzy mikrofonowej
Mf=dmax+ws;                             % długość filtru FIR
H=zeros(Mf,N);                          % macierz wag filtrów FIR
for k=1:N
    H(shf(k)+1:shf(k)+ws,k)=h;           % opóźnienia
end
%------------------------------------------------------------
for k=1:N
    y1(:,k)=filter(H(:,k),1,X(:,k));    % filtracja kanałów
end
s=ones(1,N)/N;                          % przygotowanie wektora sumującego
y=y1*s';                                % suma filtrowanych kanałów / N

Przykład 17.5

W przykładzie 17.5 (plik gcc_phat.m) przedstawiono implementację w języku MATLAB funkcji GCC-PHAT.
% Przykład 17.5
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

function y=gcc_phat(x1, x2)
%% y=gcc_phat(x1, x2)
%
% x1, x2 - sygnały wejściowe
% y - wyjście GCC-PHAT
NFFT = 2*(size(x1,1))-1;                % obliczenie długości FFT
G1=fft(x1,NFFT);                        % widmo X1(f)
G2=fft(x2,NFFT);                        % widmo X2(f)
G12=G1.*conj(G2);                       % iloczyn X1(f)X2(f)*
den=max(abs(G12),1e-6);                 % moduł |X1(f)X2(f)*|
G=G12./den;                             % X1(f)X2(f)*/|X1(f)X2(f)*|
g=real(ifft(G));                        % transformacja w dziedzinę czasu
y=fftshift(g);

Przykład 17.6

W przykładzie 17.6 (plik tdoa.m) przedstawiono implementację w języku MATLAB funkcji realizującej metodę TDOA dla kilku źródeł dźwięku.
% Przykład 17.6
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

function theta=tdoa(x1, x2, fp, l, M)
%% theta=tdoa(x1, x2, fp, l, M)
%
% x1, x2 - sygnały wejściowe, fp - częstotliwość próbkowania
% l - odległość mikrofonów, M - liczba źródeł dźwięku
v=340;                                  % prędkość dźwięku [m/s]
xx=gcc_phat(x1,x2);                     % korelacja GCC-PHAT
[val,i]=sort(xx,'descend');             % sortowanie M-maksimów
S=size(x1,1);                           % długość sygnału
for k=1:M                               % iteracja po źródłach dźwięku
    shf=S-i(k);                         % obliczenie kąta wiązki
    d=shf/fp;
    alfa=asin((d*v)/(l));
    theta(k)=alfa*180/pi;
end

Przykład 17.7

W przykładzie 17.7 (plik przyklad_17_7.m) przedstawiono program w języku MATLAB, który umożliwia śledzenie źródła dźwięku z wykorzystaniem metody TDOA.
% Przykład 17.7
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

clear all; clc;
l=0.08;                                 % odległość mikrofonów w macierzy
k=3;                                    % kanał do obliczenia TDOA względem kanału 1
[X, fp] = wavread('mch_truck.wav');     % odczyt pliku
NFFT=16384;                             % rozmiar okna FFT
K=size(X,1);                            % liczba próbek w pojedynczym kanale
wn = window(@hann,NFFT);                % okno dla FFT
j=1;
for(i=0:(NFFT/2):K-NFFT)                % przesuwanie okna z zakładką 1/2 długości okna
    XW(1:NFFT,:)=X(i+1:i+NFFT,:);       % wycinek sygnału o długości okna
    xw1=wn.*XW(:,1);                    % przemnożenie 1 kanału przez okno (wn)
    xw2=wn.*XW(:,1+k);                  % przemnożenie k-tego kanału przez okno (wn)
    theta(j)=tdoa(xw1, xw2, fp, k*l, 1);% obliczenie TDOA dla 1 źródła sygnału
    j=j+1;
end
tm=linspace(0,j*(NFFT/2)/fp,j-1);       % przestrzeń czasu
plot(tm, theta,'b*-');
xlabel(['Czas (s)']); ylabel(['Kąt (stopnie)']); grid on;

Przykład 17.8

W przykładzie 17.8 (plik przyklad_17_8.m) przedstawiono implementację w języku MATLAB wykorzystania algorytmu Filter-Sum do poprawy SNR.
% Przykład 17.8
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

clear all; clc;                             % przygotowanie obliczeń
%------------------------ Parametry macierzy ------------------------------
l=0.02;                                     % odległość mikrofonów w macierzy
%theta=30;                                  % kąt wiązki (-90 - +90)
%--------------------------------------------------------------------------
[X, fp] = wavread('mch_noise.wav');
%------------------------------ TDOA --------------------------------------
N=size(X,2);                               % liczba kanałów
theta=tdoa(X(:,1), X(:,N), fp, (N-1)*l, 1) % obliczenie kąta na podstawie TDOA
%---------------------- Narrowband Beamforming ----------------------------
%wagi filtru pasmowo-przepustowego (Hamming 10rząd 2kHz-16kHz)
hnb=[-0.00965926 0 -0.0309026 -0.153538 0.18252 0.603916 0.18252 -0.153538 -0.0309026 0 -0.00965926];
%----------------------- Broadband Beamforming ----------------------------
%wagi dla delay-sum [1 0 0 0 ...]
hbb=[1];
%--------------------------------------------------------------------------
fsum=filter_sum(X, fp, hnb, l, theta);      % filter-sum
%fsum=delay_sum(X, fp, l, theta);           % delay-sum
%--------------------------- Spektrogram ----------------------------------
amin=1000;                                  % dynamika spektrogramu

[C,f,t]=specgram(X(:,1),1024,fp,256,192);   % pojedynczy kanał
subplot(2,1,1);
bmin=max(max(abs(C)))/amin;
imagesc(t,f,20*log10(max(abs(C),bmin)/bmin));
axis xy; xlabel('Czas (s)'); ylabel('Częstotliwość (Hz)'); colorbar;

[C,f,t]=specgram(fsum,1024,fp,256,192);     % wyjście filter-sum
subplot(2,1,2);
bmin=max(max(abs(C)))/amin;
imagesc(t,f,20*log10(max(abs(C),bmin)/bmin));
axis xy; xlabel('Czas (s)'); ylabel('Częstotliwość (Hz)'); colorbar;
%--------------------------------------------------------------------------
wavwrite(fsum, fp, 'fs_noise.wav');         % zapis wyjścia do pliku
sound(fsum, fp);                            % odtwarzanie wyjścia

Przykład 17.9

W przykładzie 17.9 (plik przyklad_17_9.m) przedstawiono implementację w języku MATLAB wykorzystania algorytmu Filter-Sum do separacji źródeł dźwięku.
% Przykład 17.9
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

clear all; clc;                             % przygotowanie obliczeń
%------------------------ Parametry macierzy ------------------------------
l=0.02;                                     % odległość mikrofonów w macierzy
%thetaA=30;                                 % kąt wiązki A (-90 - +90)
%thetaB=0;                                  % kąt wiązki B (-90 - +90)
%thetaC=-45;                                % kąt wiązki C (-90 - +90)
%--------------------------------------------------------------------------
[X, fp] = wavread('mch_separation.wav');
%------------------------------ TDOA --------------------------------------
N=size(X,2);                                % liczba kanałów w pliku
theta=tdoa(X(:,1), X(:,N), fp, (N-1)*l, 3); % obliczenie 3 kątów na podstawie TDOA
thetaA=theta(1)
thetaB=theta(2)
thetaC=theta(3)
%----------------------- Narrowband Beamforming ---------------------------
%wagi filtru pasmowo-przepustowego (Hamming 10rzad 2kHz-16kHz)
hnb=[-0.00965926 0 -0.0309026 -0.153538 0.18252 0.603916 0.18252 -0.153538 -0.0309026 0 -0.00965926];
%----------------------- Broadband Beamforming ----------------------------
%wagi dla delay-sum [1 0 0 0 ...]
hbb=[1];
%--------------------------------------------------------------------------
fsum1=filter_sum(X, fp, hnb, l, thetaA);
fsum2=filter_sum(X, fp, hnb, l, thetaB);
fsum3=filter_sum(X, fp, hnb, l, thetaC);
%--------------------------- Spektrogram ----------------------------------
amin=1000;                                  % dynamika spektrogramu

[C,f,t]=specgram(X(:,1),1024,fp,256,192);   % pojedynczy kanał
subplot(2,2,1);
bmin=max(max(abs(C)))/amin;
imagesc(t,f,20*log10(max(abs(C),bmin)/bmin));
axis xy; xlabel('Czas (s)'); ylabel('Częstotliwość (Hz)'); colorbar;

[C,f,t]=specgram(fsum1,1024,fp,256,192);     % wyjście filter-sum 1
subplot(2,2,2);
bmin=max(max(abs(C)))/amin;
imagesc(t,f,20*log10(max(abs(C),bmin)/bmin));
axis xy; xlabel('Czas (s)'); ylabel('Częstotliwość (Hz)'); colorbar;

[C,f,t]=specgram(fsum2,1024,fp,256,192);     % wyjście filter-sum 2
subplot(2,2,3);
bmin=max(max(abs(C)))/amin;
imagesc(t,f,20*log10(max(abs(C),bmin)/bmin));
axis xy; xlabel('Czas (s)'); ylabel('Częstotliwość (Hz)'); colorbar;

[C,f,t]=specgram(fsum3,1024,fp,256,192);     % wyjście filter-sum 3
subplot(2,2,4);
bmin=max(max(abs(C)))/amin;
imagesc(t,f,20*log10(max(abs(C),bmin)/bmin));
axis xy; xlabel('Czas (s)'); ylabel('Częstotliwość (Hz)'); colorbar;
%--------------------------------------------------------------------------
wavwrite(fsum1, fp, 'fs_1.wav');            % zapis wyjscia 1 do pliku
wavwrite(fsum2, fp, 'fs_2.wav');            % zapis wyjscia 2 do pliku
wavwrite(fsum3, fp, 'fs_3.wav');            % zapis wyjscia 3 do pliku
sound(fsum1, fp);                            % odtwarzanie wyjscia 1
sound(fsum2, fp);                            % odtwarzanie wyjscia 2
sound(fsum3, fp);                            % odtwarzanie wyjscia 3

Przykład 17.10

W przykładzie 17.10 (plik delay_speaker.m) przedstawiono implementację w języku MATLAB funkcji sterującej wiązką w macierzy głośnikowej.
% Przykład 17.10
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

function Y=delay_speaker(x, fp, N, l, theta)
%% Y=delay_speaker(x, fp, N, l, theta)
%
% x - sygnał wejściowy, fp - częstotliwość próbkowania
% N - liczba kanałów, l - odległość głośników, theta - kąt wiązki
% Y - wielokanałowy sygnał wyjściowy
v=340;                                  % prędkość dźwięku
alfa=(theta*pi)/180;                    % stopnie/radiany
d=(l*sin(alfa))/v;                      % opóźnienie [s]
M=size(x,1);                            % liczba próbek w kanale
for k=1:N
    shf(k)=floor((k-1)*abs(d)*fp);      % wektor opóźnień dla
end                                     % głośników
if(d<0)                                 % transformacja dla
    shf=fliplr(shf);                    % ujemnego kąta wiązki
end
Y = zeros(M,N);
for k=1:N                               % wiązka sterowana
    Y(1+shf(k):M,k)=x(1:M-shf(k));
end

Przykład 17.11

W przykładzie 17.11 (plik przyklad_17_11.m) przedstawiono implementację w języku MATLAB algorytmu dla macierzy głośnikowej z trzema niezależnymi wiązkami sterowanymi.
% Przykład 17.11
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 17. Macierze mikrofonowe i głośnikowe. 
% Autor rozdziału: Daniel Król.

clear; clc;                             % przygotowanie obliczeń
%------------------------ Parametry macierzy ------------------------------
l=0.02;                                 % odległość głośników
N=16;                                   % liczba kanałów

thetaA=30;                              % kąt wiązki A (-90 - +90)
thetaB=0;                               % kąt wiązki B (-90 - +90)
thetaC=-45;                             % kąt wiązki C (-90 - +90)
%------------------------ Sygnały wejściowe -------------------------------
[x0, fp] = wavread('bird1.wav');        % sygnał A
[x1, fp] = wavread('bird2.wav');        % sygnał B
[x2, fp] = wavread('bird3.wav');        % sygnał C
%----------------------- Generowanie wiązek -------------------------------
Y0=delay_speaker(x0, fp, N, l, thetaA);
Y1=delay_speaker(x1, fp, N, l, thetaB);
Y2=delay_speaker(x2, fp, N, l, thetaC);
%--------------------------------------------------------------------------
Y=(Y0+Y1+Y2)/3;                         % znormalizowana suma wiązek
wavwrite(Y, fp, 'mch_separation.wav');  % zapis wyjścia