function synthesized_signal = Spectral_Splice_Repeat(magnitude_signal,phase_signal,window_size,overlap)
%%
% This function creates a signal with the magnitude spectrum of one input
% signal and the phase spectrum of another. 
%
% spectral_splice_repeat breaks each input into segments whose length is
% defined by the input 'window_size. It transforms one segment from each
% signal one-at-a-time into the time domain via the complex discrete Fourier
% transform by using Matlab's built-in fft function, which uses the equation:
%
%                    N
%      X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                   n=1
% 
% The absolute value of one segment is used to obtain its magnitude spectrum,
% the angle of the other segment is used to obtain it's phase spectrum. The
% two values are used to create a complex segment which is then transformed
% back into time domain to create a synthesised segment.  All the
% synthesised segments are added together to create a synthesized_signal
% that is equal in length to the input signals. Each segment is multiplied
% by a hanning window and overlaps the previous signal by 90 percent to
% create a smoothly changing signal. 
%
% The function spectral_splice_repeat varies from the function
% spectral_splice in that is compensates for quiet sections of the
% magnitude_signal by comparing them to a set rms threshold. If a segment of
% the magnitude_signal drops below this threshold it is discarded. In its
% place, the most recent segment of the magnitude_signal to exceed the
% threshold is used until another segment of the magnitude_signal exceeds
% the threshold. If the magnitude_signal has yet to exceed the
% threshold a sine wave generated by the function is used in its place.
%
% Inputs are: 
% magnitude_signal and phase_signal are mono audio signals whose spectral
% data will be combined.
%
% window_size - The length of the segments the signals will be broken into.
%               AKA the amount of samples to be analysed at a time
% overlap     - The amount each window overlaps the preceding one, expressed as a percentage
 
%% A segment of a sine wave is created to be used in the case that the magnitude_signal has yet to exceed the threshold
sinewave_length=0:window_size; %Creates a vector of increasing integers equal in length to the window size. 
sinewave_segment=0.1.*sin(sinewave_length./20); %Creates a sine wave equal in length to the window size.
last_fft_magnitude_signal_segment=fft(transpose(sinewave_segment)); % Ensures that sinewave_segment will be used until a magnitude_segment exceeds the rms threshold. This variable will be redefined every time a magnitude_segment exceeds the rms threshold.
 
%% The segments and the synthesized_signal are defined.
num_windows = ceil(length(magnitude_signal) / (window_size*(1-(overlap/100)))); % The number of windows
% are calculated taking the amount of overlap into account. The ceil
% function is used to round values up to the nearest integer.
wavelength = length(magnitude_signal); %Defines the length of the magnitude_signal for further calculations
synthesized_signal = zeros(wavelength, 1); % Creates a variable equal in length to the magnitude_signal (which is equal in length to the phase_signal.
% Processed signal segments will later be added to this variable.
 
%% One segment from signal is transformed into the time domain and the two segments are synthesised. This process repeats until the entire signal is synthesised.
for i = 0:num_windows-1
 window_start = (i * window_size + 1)-(window_size*i*(overlap/100)); % The beginning of the current segment is calculated taking overlap into account.
 window_end = window_start + window_size; % The end of the current segment is calculated using the size of the window_size parameter.
  if window_end > wavelength;  % This prevents any segment extending past the end of the synthesized_signal.
      window_end = wavelength;
  end
    
 current_fft_magnitude_signal_segment = fft(magnitude_signal(window_start:window_end).*hann((window_end-window_start)+1)); % Multiplies the current segment of the magnitude_signal by a hanning window to smooth its corners and  creates a frequency domain signal for the current segment of the magnitude_signal by using the fft function. 
 current_fft_phase_signal_segment = fft(phase_signal(window_start:window_end).*hann((window_end-window_start)+1)); % Multiplies the current segment of the phase_signal by a hanning window to smooth its corners and Creates a frequency domain signal for the current segment of the phase_signal by using the fft function. 
 
 r=rms(magnitude_signal(window_start:window_end)); %Finds the rms value of the current segment from the magnitude signal in order to avoid using quiet portions of the signal.
 same_size=(length(current_fft_magnitude_signal_segment)); % This length is used to prevent errors when the last_fft_magnitude_signal_segment is used
 
 %unused_phasespectrum = angle(current_fft_magnitude_signal_segment); %These values are just used for plotting and illustration
 %unused_magspectrum = abs(current_fft_phase_signal_segment); %These values are just used for plotting and illustration
 
  if r>0.04; % Sets an rms threshold for the current frequency domain magnitude_signal_segment, if the current_fft_magnitude_signal_segment is below this threshold, the function will use the last segment whose RMS value exceeded this threshold.
 
  magspectrum = abs(current_fft_magnitude_signal_segment);% % Obtain magnitude spectrum of the magnitude_signal's current segment by taking the its absolute value and disregarding any imaginary numbers. 
  last_fft_magnitude_signal_segment=current_fft_magnitude_signal_segment; % Saves the current magnitude segment as last_fft_magnitude_signal_segment every time it exceeds the set threshold
  phasespectrum = angle(current_fft_phase_signal_segment); % Obtain phase spectrum of the phase_signal's current segment by calculating its angle from the origin (theta) where phasespectrum is a series of complex numbers and theta=arctan(imaginary/real).
    
  else
  magspectrum = abs(last_fft_magnitude_signal_segment(1:same_size)); % Uses the last_fft_magnitude_signal_segment to meet the RMS threshold if the current_fft_magnitude_signal_segment is too quiet. It obtains the segment's magnitude spectrum by taking its absolute value.
  phasespectrum = angle(current_fft_phase_signal_segment); %Obtain phase spectrum of the phase_signal's current segment.
  end
    
 complexspectrum = magspectrum .* exp(1i .* phasespectrum); %The complex frequency domain segment is created using the formula complex signal = magnitude * e^(i(phase)) and the variables magspectrum and phasespectrum
 processed_audio_signal_segment = real(ifft(complexspectrum)); %Transforms the new frequency domain signal into the time domain. 
 synthesized_signal(window_start:window_end) = synthesized_signal(window_start:window_end) + processed_audio_signal_segment; %Each processed_audio_signal_segment is added to the correct time index by referencing the window_start and window_end variables.
end
  synthesized_signal=0.99.*synthesized_signal./max(abs(synthesized_signal));%Normalises the output signal.

