Перейти к содержанию
    

Beamformer(лучеформирователь) в матлаб

Имеется код Delay and Sum бимформера для двух сенсоров. Код дан как приложение к статье http://www.google.com/url?sa=t&rct=j&a...kA4pTsNACkZq0ZA

%% INITIALIZE
clear all; close all; clc;

%Properties of the phone
d = 0.12;					   % Distance between the mics
N = 2;						  % Amount of mics

% Physicalconstants
c = 343;						% Speed of sound in air

% Properties of the signal
freq = 300;					 % The frequency of the source signal
freq2 = 600;					% The frequency of the noise signal
theta = pi;					 % Angle of the source signal
theta2 = 2*pi/9;				% Angle of the noise signal

% Our chosen properties
fs = 80000;					 % Sample frequency
time_frame = 0.02;			  % Duration of 1 time frame
sptf = fs * time_frame;		 % The amount of samples per time frame
noise_amp = 0.1;				% The noise amplitude
L = 2^ (ceil(log2(sptf)));	  % The power of 2 for the FFT
totaltime = 1;				  % The total time of the signal

%% INPUT SIGNAL
tt = linspace (0, totaltime, fs * totaltime).';	   % Samples till total time
tau1 = d/c * sin(theta);							% Time delay to mic 2 for signal
tau2 = d/c * sin(theta2);						   % Time delay to mic 2 for noise
signal = sin (2 * pi * freq * tt);				  % Signal mic 1
signal2 = sin (2 * pi * freq * (tt - tau1));		% Signal mic 2

noise = 4 * sin(2 * pi * freq2 * tt) + noise_amp * randn(totaltime * fs, 1);			% Noise mic 1
noise2 = 4 * sin(2 * pi * freq2 * (tt - tau2)) + noise_amp * randn(totaltime * fs, 1);  % Noise mic 2

signal_total = signal + noise;	  % Sum of the signal and noise at mic 1
signal_total2 = signal2 + noise2;   % Sum of the signal and noise at mic 2
x1 = [signal, signal2];
xsum = [signal_total, signal_total2];

f_center = linspace(0, fs - fs/(2 * L), L).';		 % The center frequencies for each bin
zeta = -1i * 2 * pi * f_center * d * sin(theta)/c;	  % The phase shift of the signal for mic 2
a_n = 1/N;											  % The amplitude weights
w = a_n * [(exp(zeta)), ones(L, 1)];					% Weights with delays

%% RECONSTRUCTING THE ORIGINAL SIGNAL
nr_frames = floor((length(signal)-sptf)/(sptf/2));
rec_signal = zeros(size(signal(:, 1)));
rec_signal_total = zeros(size(signal_total(:, 1)));

testsignal1 = zeros(L, nr_frames);
testsignal2 = zeros(L, nr_frames);

for I = 1 : nr_frames
win = [sqrt(hanning(sptf)), sqrt(hanning(sptf))];
frame_x1 = x1((sptf/2)*(I-1)+1:(sptf/2)*(I-1)+sptf , :).*(win);
fft_x1 = fft (frame_x1, L);

frame_xsum = xsum((sptf/2)*(I-1)+1:(sptf/2)*(I-1)+sptf, :).*(win);
fft_xsum = fft(frame_xsum, L);
fft_xsum (end, :) = real(fft_xsum(end, :));

fft_x1([1, L/2 + 1], :) = real(fft_x1([1, L/2 + 1], :));
part1 = fft_x1(1 : L/2 + 1, :);
fft_x1 = [part1; conj(flipud(part1(2 : end - 1, :)))];

fft_xsum ([1, L/2 + 1], :) = real(fft_xsum([1, L/2 + 1], :));
part1 = fft_xsum (1 : L/2 + 1, :);
fft_xsum = [part1; conj(flipud(part1 (2 : end - 1, :)))];
estimate_fft_signal = w .* (fft_x1);
estimate_fft_signal_2 = estimate_fft_signal(:, 1) + estimate_fft_signal(:, 2);

estimate_fft_signal_total = w .* (fft_xsum);
estimate_fft_signal_total_2 = estimate_fft_signal_total(:, 1) + estimate_fft_signal_total(:, 2);

testsignal2(: , I) = estimate_fft_signal_total_2;
testsignal1(: , I) = fft_xsum(: , 1);

estimate_signal = real(ifft(estimate_fft_signal_2));

estimate_signal_total = real (ifft(estimate_fft_signal_total_2));
rec_signal((sptf / 2) * (I - 1) + 1 : ( sptf / 2) * (I - 1) + sptf ) = rec_signal((sptf / 2) * (I - 1) + 1 : (sptf/2) * (I - 1) + sptf) + estimate_signal(1 : sptf) .* (sqrt(hanning(sptf)));
rec_signal_total((sptf / 2) * (I - 1) + 1 : (sptf / 2) * (I - 1) + sptf) =  rec_signal_total((sptf / 2) * (I - 1) + 1 : (sptf / 2) * (I - 1) + sptf) +  estimate_signal_total(1 : sptf) .* (sqrt(hanning(sptf)));

end

%% PLOTS
% Time plot
t = linspace(0 , totaltime , totaltime * fs);
figure;
subplot (2, 1, 2)
plot (t, real(signal_total2(1 : totaltime * fs)));
hold on
plot (t , real( rec_signal_total(1 : totaltime * fs)), '.-r');
plot (t , real( signal2(1 : totaltime * fs)), '-g');
legend ('nois y in', 'beamformed', 'desiredin')
title ('Incoming, beamformed and desired signals')
axis ([0.2000 0.210 -5 5]);
xlabel('Time(inseconds)');
ylabel('Amplitude');

% Fr equency plot
subplot (2, 1, 1)
plot (linspace(1, fs, L), 10 * log10(mean(abs(testsignal1).^2, 2)));
hold on
plot (linspace(1, fs, L), 10 * log10 (mean(abs(testsignal2).^2, 2)), 'r');
title ('Spect rum plot s of in signal and signal after beamforming i n dB')
legend ('insignal', 'beamformed signal')
axis ([0 4000 10 50]);
xlabel ('Fr equency(inHertz)');
ylabel ('Magnitude(in dB)');

% Polarplot
figure;

polarfreq = [300, 600];
polarbin = zeros(1, length(polarfreq));

for I = 1 : length (polarfreq);
polarbin (I) = floor(polarfreq(I)/fs * L);

delta = ((c / f_center(polarbin(I))) / d )^ -1; 
[w_dakje, polarhoek] = beam_resp(w(polarbin(I), :), L, delta);

subplot (ceil(length(polarfreq) / 2), 2, I)
polar90 (polarhoek, abs(w_dakje));
title (['Fr equency =', num2str (polarfreq(I))])

end

Матлаб на % Polarplot выдает ошибку типа

Undefined function 'beam_resp' for input arguments of type 'double'.

Error in DelayAndSum (line 122)
    [w_dakje, polarhoek] = beam_resp(w(polarbin(I), :), L, delta);

Подскажите пожалуйста, что не так и как исправить.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Функция не определена. Видимо, нужно просить автора статьи предоставить код этой функции.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

Вы можете написать сейчас и зарегистрироваться позже. Если у вас есть аккаунт, авторизуйтесь, чтобы опубликовать от имени своего аккаунта.

Гость
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
  • Создать...