clear
clc

pause on


%%  Read in the sound data (d = data, r = sample rate)
%[d,r] = audioread('LoveGameLadyGaga.ogg.wav');
%[d,r] = audioread('twinkle.wav');
%[d,r] = audioread('thank u, next short.wav');
[d,r] = audioread('Cut To The Feeling.wav');
%[d,r] = audioread('02 Blank Space.m4a');
%[d,r] = audioread('Ill Never Love Again.wav');
d = d(:,1); % Set d equal to only its first column

%sound(d,r); % Listen to it

%% Pass it through the filter
%df = bandpass(d,[196 1047],r);
df = bandpass(d,[80 250],r);
%df = bandpass(d,[500 900],r);
%df = d;
%sound(df,r);
%% Create a new spectrogram

[B,f,t] = specgram(df,1024,r);

%% Create a novelty Curve
B = abs(B(2:7,:)); % For frequencies 50 - 150
%B = abs(B(33:62,:)); % For frequencies 500 - 900

% Calculate the difference across rows of the specgram
B_diff = diff(B,1,2);

% Calculate the accumulated difference in each column
Acc = sum(B_diff,1);
Acc = Acc./max(Acc);
Acc(Acc<0) = 0;

%% Plot the Novelty Curve vs the Spectrogram
figure();
subplot(2,1,1);
specgram(df,1024,r);
title('Spectrogram');
subplot(2,1,2);
plot(t(1:length(Acc),1),abs(Acc));
axis([0 25 0 1]);
title('Novelty Curve');

%% Create a figure with 4 Novelty Curves w different averaging windows
figure();
subplot(4,1,1);
Acc = sum(B_diff,1);
Acc = Acc./max(Acc);
Acc(Acc<0) = 0;
plot(t(1:length(Acc),1),abs(Acc));
title('Novelty Curve');

Moving_Average = movmean(Acc,[5 5]);
Acc_Denoised = Acc - Moving_Average;
Acc_Denoised(Acc_Denoised<0) = 0;
subplot(4,1,2);
plot(t(1:length(Acc),1),abs(Acc_Denoised));
title('Novelty Curve Denoised (W = 11)');

Moving_Average = movmean(Acc,[10 10]);
Acc_Denoised = Acc - Moving_Average;
Acc_Denoised(Acc_Denoised<0) = 0;
subplot(4,1,3);
plot(t(1:length(Acc),1),abs(Acc_Denoised));
title('Novelty Curve Denoised (W = 21)');

Moving_Average = movmean(Acc,[25 25]);
Acc_Denoised = Acc - Moving_Average;
Acc_Denoised(Acc_Denoised<0) = 0;
subplot(4,1,4);
plot(t(1:length(Acc),1),abs(Acc_Denoised));
title('Novelty Curve Denoised (W = 51)');

%% Create a figure with 4 Novelty Curves w different averaging windows
figure();
subplot(4,1,1);
Acc = sum(B_diff,1);
Acc = Acc./max(Acc);
Acc(Acc<0) = 0;
plot(t(1:length(Acc),1),abs(Acc));
title('Novelty Curve');

subplot(4,1,2);
Acc(Acc<0.125) = 0;
plot(t(1:length(Acc),1),abs(Acc));
title('Novelty Curve Denoised (Acc<0.125 = 0)');

subplot(4,1,3);
Acc(Acc<0.25) = 0;
plot(t(1:length(Acc),1),abs(Acc));
title('Novelty Curve Denoised (Acc<0.25 = 0)');

subplot(4,1,4);
Acc(Acc<0.375) = 0;
plot(t(1:length(Acc),1),abs(Acc));
title('Novelty Curve Denoised (Acc<0.375 = 0)');

%% Plot the Novelty Curve with (w = 51 & Acc<150 = 0)
Acc = sum(B_diff,1);
Acc = Acc./max(Acc);
Acc(Acc<0) = 0;

figure();
subplot(4,1,1);
plot(t(1:length(Acc),1),abs(Acc));
title('Novelty Curve');

subplot(4,1,2);
Moving_Average = movmean(Acc,[25 25]);
Acc_Denoised = Acc - Moving_Average;
Acc_Denoised(Acc_Denoised<0) = 0;
plot(t(1:length(Acc),1),abs(Acc_Denoised));
title('Novelty Curve Denoised (W = 51)');

subplot(4,1,3);
Acc(Acc<0.30) = 0;
plot(t(1:length(Acc),1),abs(Acc));
title('Novelty Curve Denoised (Acc<0.30 = 0)');

subplot(4,1,4);
Acc_Denoised(Acc_Denoised<0.30) = 0;
plot(t(1:length(Acc),1),abs(Acc_Denoised));
title('Novelty Curve Denoised (w = 51 & Acc<0.30 = 0)');


%% Beat Plotter

Acc = sum(B_diff,1);
Acc = Acc./max(Acc);
Acc(Acc<0) = 0;

Acc_Denoised(Acc_Denoised>0) = 1;
figure()
plot(Acc_Denoised);
title('Beat Plot');
xlabel('Samples');
ylabel('Normalized Amplitude');

% Convolve the Signal with an Exponential
t1 = 0:0.2:20;
decay = exp(-t1./4);
intensity = conv(Acc_Denoised,decay);
intensity(intensity>1) = 1;
intensity = fix(255.*intensity);

% Convolve the Signal with a Sine
%t1 = -pi:0.1:pi;
% sine = (cos(t1)+1)./2;
% intensity = conv(Acc_Denoised,sine);
% intensity(intensity>1) = 1;

rgb = zeros(length(intensity),3);
rgb(:,1) = intensity;
rgb(:,2) = 0;
rgb(:,3) = 0;

%% FOR DEBUGGING
% figure();
% time = linspace(0,length(intensity),length(intensity));
% plot(time,intensity/255);
% title('Beat Plot Convolved with Partial Sine Wave');
% xlabel('Samples');
% ylabel('Normalized Amplitude');

%% Plot the data real-time

% figure();
% h = plot(t(1:length(Acc),1),abs(Acc_Denoised));
% h.YData = zeros(1,length(Acc_Denoised));
% axis([0 24 0 500]);
% y_temp = zeros(1,length(Acc_Denoised));
% for n = 1:length(Acc_Denoised)
%     y_temp(1,n) = Acc_Denoised(1,n);
%     h.YData = y_temp;
%     drawnow
%     %pause(0.01);
% end