% This is the main code. 

clc 
clear
close all
format long
set(0,'DefaultTextInterpreter','none'); % this is so underscores in file names aren't interpreted as subscripts

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%    CUSTOMIZATION SECTION     %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% ENTER FID FILENAME HERE %%%
filename = "Glutamate_FID.csv";

FID = readtable(filename); 
FIDy = FID{:, 2};
maxSamples = length(FIDy);
samplesX = [1:maxSamples];
maxf = 2050;                % frequencies will be resynthesized from 0 - "maxf" Hz. 
Fs = 4096;                  % sampling frequency
duration = maxSamples/Fs;   % duration of sound in seconds
chunknumber = 5000;         % resolution of amplitude envelope for reshaping; higher = more accurate

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fourierY = abs(fft(FIDy));

% now we must remove the second half of the fourier transform, which is a
% redundant mirror image.
killmat = [1:(maxSamples/2)] + (maxSamples/2);
fourierY(killmat) = [];

% to maintain the number of samples, we must interpolate the data to
% regenerate the matrix of size "maxSamples". 

queryX = 0.5 : 0.5 : maxSamples/2;
fourierY = interp1([1:length(killmat)], fourierY, queryX);

normy = fourierY/max(fourierY);

% modify normy to preserve the second and penultimate-index values
% (findpeaks function won't see the first or last values)
% through the peak-picking process by setting to 1, so they will be kept.
% This is so that the precise scale and location of the frequency axis is
% maintained. 

normy(1) = 0;
normy(2) = 1;
normy(maxSamples-1) = 1;

[peakamps, peakfreqs] = findpeaks(normy);

peakfreqs = flipud(peakfreqs);

% removing all insignificant peaks (below 10% of max amplitude)
indices = find(peakamps<0.1);
peakamps(indices)=[];
peakfreqs(indices)=[];

% now we have to remove the border peaks and scale the x-axis. 
peakamps(1) = 0;
peakamps(length(peakamps)) = 0;
scaledpeakfreqs = (peakfreqs/max(peakfreqs)).*maxf;

sinlength = [(1/Fs):1/Fs:duration];
buildwave = sin(2*pi*0.*sinlength);

% This is the additive synthesis; a loop that builds a waveform out of
% sinewaves of amplitude and frequency detailed by the picked peaks

for index = 1:length(peakamps)
    buildwave = buildwave + peakamps(index)*sin(2*pi*(scaledpeakfreqs(index))*sinlength);
end

% to make an envelope for the waveform,
% analyse the waveform in "chunknumber" chunks, where each chunk is [#ofSamples / chunknumber].
% extract the max value from each chunk, 
% put it into a matrix of length n.

% First, we must ensure the file has enough samples to evenly divide by chunk size.
% If the number of samples is not divisible by the number of desired
% chunks, the code will append enough 0's to the end of the data to make a
% divisible number of samples. 
% This allows any integer size of chunk to be used.

chunkSize = ceil(maxSamples/chunknumber);

if mod(maxSamples, chunknumber)~= 0
    FIDy(chunkSize*chunknumber) = 0;
end 

chunkd = reshape(FIDy, [chunkSize, chunknumber]);

chunkMaxes = max(abs(chunkd));
chunkX = [1:chunknumber].*chunkSize;
env = interp1(chunkX, chunkMaxes, samplesX);

resynth = env.*buildwave;
resynth = resynth/(max(abs(resynth)));

figure
subplot(2,1,1);
plot(FIDy/max(FIDy));
title('original waveform of ' + filename);
xlabel('sample number');
ylabel('amplitude');
subplot(2,1,2);
plot(resynth);
xlabel('sample number');
ylabel('amplitude');
title('resynth of ' + filename)

figure
stem(scaledpeakfreqs, peakamps);
title('cleaned up Frequency Spectrum of ' + filename)
xlabel('frequency (Hz)');
ylabel('amplitude');

soundsc(resynth, Fs);
% audiowrite('output.wav', resynth, Fs);

%% Just the molecule sound

clc 
clear
format long

% EDIT FILENAME HERE
filename = "glutamate_FID.csv";

FID = readtable(filename);

Fs = 4096;
FIDy = FID{:, 2};

soundsc(FIDy, Fs);

FIDy = FIDy/(max(abs(FIDy)));
% audiowrite('moleculeSound.wav', FIDy, Fs);





