%% 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

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);