% Section 3.2(a)
% Read the contents of the audio file
waveFilename = 'violin-C4.wav';
[instrumentSound, fs] = audioread(waveFilename);

% Play back the recording with automatic scaling
soundsc(instrumentSound, fs);
 
% Plot the waveform in the time domain
N = length(instrumentSound);
Ts = 1/fs;
Tmax = (N-1)*Ts;
t = 0 : Ts : Tmax;
figure;
plot(t, instrumentSound);
xlabel('Time [s]');
ylabel('Signal amplitude');

% Section 3.2(b)
mean(instrumentSound)

% Section 3.2(c)
xlim( [1.5 1.6] );

fundPeriod = (1.599 - 1.50295) / 25
fundFrequency = 1 / fundPeriod
C4fundFrequency = 261.63
fundFrequencyRelErr = (C4fundFrequency - fundFrequency) / C4fundFrequency

% Section 3.3(a)
% Plot the magnitude of the frequency content
% using the discrete-time Fourier series
fourierSeriesCoeffs = fft(instrumentSound);
N = length(instrumentSound);
freqResolution = fs / N;
ff = (-fs/2) : freqResolution : (fs/2)-freqResolution;
figure;
plot(ff, abs(fftshift(fourierSeriesCoeffs)));
xlabel('f');

% Section 3.3(b)
N = length(instrumentSound);
nonNegativeFrequency = 261.45;
f0 = fs / N;
fftIndex = round(nonNegativeFrequency / f0) + 1;   %%% Value is 914

abs(fourierSeriesCoeffs(fftIndex))
angle(fourierSeriesCoeffs(fftIndex))

[maxval, maxind] = max(abs(fourierSeriesCoeffs));
gainValue = abs(fourierSeriesCoeffs(maxind))
phaseValue = angle(fourierSeriesCoeffs(maxind))
freqValue = (maxind - 1)*fs/N


% Section 3.3(c)
% Plot the spectrogram
figure;
blockSize = round(N/4);
overlap = round(0.875 * blockSize);
spectrogram(instrumentSound, blockSize, overlap, blockSize, fs, 'yaxis');

powVal = spectrogram(instrumentSound, blockSize, overlap, blockSize, fs, 'yaxis');
[maxValue, maxIndex] = max(abs(powVal), [], 'all', 'linear');
[row, col] = ind2sub(size(powVal), maxIndex);
maxPowVal = abs(powVal(row, col));

% Section 4.0(a)
synthFourierSeriesCoeffs = zeros(N, 1);

synthFourierSeriesCoeffs(914) = fourierSeriesCoeffs(914);
synthFourierSeriesCoeffs(37588) = fourierSeriesCoeffs(37588);
synthTwoPeaks = ifft(synthFourierSeriesCoeffs);
soundsc(synthTwoPeaks, fs);
pause(4);

synthFourierSeriesCoeffs(1827) = fourierSeriesCoeffs(1827);
synthFourierSeriesCoeffs(36675) = fourierSeriesCoeffs(36675);
synthFourPeaks = ifft(synthFourierSeriesCoeffs);
soundsc(synthFourPeaks, fs);
pause(4);

% Section 4.0(b)
% Needs fourierSeriesCoeffs vector computed in Section 2.3 above
N = length(fourierSeriesCoeffs);
fourierSeriesCoeffsAbs = abs(fourierSeriesCoeffs);

% Nkeep must be even to have an equal number of negative and positive freq. 
Nkeep = 10;
frequenciesKept = zeros(Nkeep, 1);

% Find the Nkeep strongest positive and negative frequency components
synthSoundCoeffs = zeros(N, 1);
for n = 1:Nkeep
    [ai, i] = max(fourierSeriesCoeffsAbs);
    synthSoundCoeffs(i) = fourierSeriesCoeffs(i);
    fourierSeriesCoeffsAbs(i) = 0;

    k = i - 1;
    if ( k < N/2 )
       frequenciesKept(n) = fs*k/N;
    else
	   frequenciesKept(n) = fs*((k/N) - 1);
    end
end

% Convert Fourier series coefficients to time domain using inverse FFT
synthSound = ifft(synthSoundCoeffs);    
soundsc(synthSound, fs);
frequenciesKept

% Section 4.0(c)
% Needs fourierSeriesCoeffs vector computed in Section 2.3 above
N = length(fourierSeriesCoeffs);
fourierSeriesCoeffsAbs = abs(fourierSeriesCoeffs);

% Nkeep must be even to have an equal number of negative and positive freq. 
Nkeep = 10;
frequenciesKept = zeros(Nkeep, 1);

% Frequency masking +/- 7.5 Hz around a center frequency
freqOffset = 7.5;
indexOffset = round(freqOffset*N/fs);

% Find the Nkeep strongest positive and negative frequency components
synthSoundCoeffs = zeros(N, 1);
for n = 1:Nkeep
    [ai, i] = max(fourierSeriesCoeffsAbs);
    synthSoundCoeffs(i) = fourierSeriesCoeffs(i);

    % Frequency masking +/- 7.5 Hz around a center frequency
    % Keep indexing in fourierSeriesCoeffsAbs in range of [1, N]
    imin = max(i - indexOffset, 1);
    imax = min(i + indexOffset, N);
    fourierSeriesCoeffsAbs(imin:imax) = 0;

    k = i - 1;
    if ( k < N/2 )
       frequenciesKept(n) = fs*k/N;
    else
	   frequenciesKept(n) = fs*((k/N) - 1);
    end
end

% Convert Fourier series coefficients to time domain using inverse FFT
synthSound = ifft(synthSoundCoeffs);    
soundsc(synthSound, fs);
frequenciesKept
