%% Load wave file 
% Wave files store additional information about
% sampling rate, usually left and right channel amplitudes, how many bits per sample etc. 
% Matlab knows how to extract it.
[data, SampleRate, NbitsPerSample]=wavread('voice_record.wav');

y=data(:,1); % we will work only with one channel in this demo 
N=length(y); %number of sample points
% if you need other one use data(:,2)
t=( (1:N)*1/SampleRate ).'; % time spacing is 1/SampleRate and we want column vector

%% %%%%%%%%%%%%%%% Generally you  have y vs t data set %%%%%%%%%%%%%%%%%%%%%%%%%%%
% see how the processing of such is done in the commented block 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%___%%% Load data file
%___%data=load('y_vs_t.dat'); 
%___%
%___%t=data(:,1);
%___%y=data(:,2);
%___%N=length(t);
%___%
%___%%% Let's find out the sample rate
%___%dt_per_point=diff(t);  % notice that matlab has this shortcut dt(i)=t(i+1)-t(i)
%___%
%___%% First, look at the data points spacing
%___%figure(2);
%___%plot(dt_per_point);
%___%title('Data spacing vs point index');
%___%xlabel('Data index');
%___%ylabel('Time (S)');
%___%
%___%%% If spacing looks almost equidistant, we just approximate it with mean.
%___%   Otherwise we need to interpolate date, over equidistant points
%___%dt=mean(dt_per_point);
%___%
%___%SampleRate=1/dt;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Let's plot data
figure(1)
hold off;
niceplotset;
plot(t,y);
title('Signal vs time');
xlabel('Time (S)');
ylabel('Amplitude');

%% Let's hear the sound. Sometimes matlab does not do it right
% sound(y,SampleRate);

%% Now we calculate Fourier transform of the data
Y=fft(y);

%% Calculate corresponding frequency for the DFT
spectrum_freq=fourier_frequencies(SampleRate, N);

%% Plot the spectrum amplitudes 
figure(2); hold off;
niceplotset;
% Since we have negative frequencies at the right, we need to find index of sorted freq
[dummy,index]=sort(spectrum_freq);
plot(spectrum_freq(index), abs(Y(index)) ); % generally Y is complex, so we plot |Y|
xlabel('Frequency (Hz)');
ylabel('FFT amplitude');
title('Spectrum amplitudes');
% Have a look at this plot and chose proper the filter roll off frequency

%% Let's filter the data with a low-pass filter
Freq3db=100; % The low-pass filter 3dB point.
[spectrum_filtered,g_vs_freq]=freq_space_low_pass_filter(Y, SampleRate, Freq3db);

%% Show filter's characteristics with a Bode plot
figure(3);
subplot(2,1,1);
niceplotset;
plot( spectrum_freq(index), abs(g_vs_freq(index)) );
xlabel('Frequency (Hz)');
ylabel('Filter gain');

subplot(2,1,2);
niceplotset;
plot( spectrum_freq(index), angle(g_vs_freq(index)) );
xlabel('Frequency (Hz)');
ylabel('Filter phase');

%% Show filtered and unfiltered spectra
figure(4); hold off;
niceplotset;
plot( spectrum_freq(index), abs(Y(index)), 'b' ); hold on;
plot( spectrum_freq(index), abs(spectrum_filtered(index)), 'r' ); 
legend('raw', 'filtered');
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude');

%% Reconstruct the filtered signal
y_filtered=ifft(spectrum_filtered);

%% Check that the filter works correct, i.e., there is no significant imaginary part
figure(5); hold off;
niceplotset;
plot(t,imag(y_filtered));
title('Filtered signal imaginary part vs time');
xlabel('Time (S)');
ylabel('Amplitude');

%% Get rid of imaginary part if above justifies it.
y_filtered=real(y_filtered);

%% Compare the filtered signal and the original raw one.
figure(1); hold off;
niceplotset;
plot(t, y, 'b-'); hold on;
plot(t, y_filtered,'r-');
title('Signal vs time');
xlabel('Time (S)');
ylabel('Amplitude');
legend('raw', 'filtered');

%% Plat the filtered sound wave.
% sound(y_filtered, SampleRate);

%% Save the filtered wave into a file.
wavwrite(y_filtered/(0.1+max(y_filtered)), SampleRate, NbitsPerSample, 'voice_record_filtered.wav');


%% Down sampling.
% The high frequency part of the filtered  spectrum has very small amplitude.
% We will truncate the spectrum, and keep only frequencies with in the filter band.
% This way we reduce number of points in the spectrum and thus in the reconstructed
% time domain signal as well.
%
% There is some flexibility in a choice of the cut off frequency. 
% It depends on the filter properties. For a single stage RC or RL type filter,
% the 10*Freq3db point is a reasonable choice.
indexes=(abs(spectrum_freq) < 10*Freq3db);
spectrum_freq_down_sampled = spectrum_freq(indexes);
spectrum_down_sampled = spectrum_filtered(indexes);
N_down_sampled = length(spectrum_down_sampled);
% we need to  normalize the spectrum
spectrum_down_sampled=spectrum_down_sampled*N_down_sampled/N;

SampleRate_down_sampled=SampleRate*N_down_sampled/N;

y_down_sampled=real(ifft(spectrum_down_sampled));
t_down_sampled = ( (1:N_down_sampled)*1/SampleRate_down_sampled ).';


%% compare filtered and raw signals
figure(1); hold off;
niceplotset;
%plot(t, y, 'b-'); hold on;
%legend('raw');
plot(t, y_filtered,'r-'); hold on;
plot(t_down_sampled, y_down_sampled,'k-','linewidth',2);
title('Signal vs time');
xlabel('Time (S)');
ylabel('Amplitude');
legend('filtered', 'down sampled');
% you probably will not able to distinguish filtered and down sampled version

% Note that down sampling and decimation (i.e., keeping every Nth point)
% are not equivalent procedures. Down sampling is removing high frequency components, 
% while decimation move them to a lower frequency domain (recall the reflection phenomenon).
% Do not use the decimation unless you know what you are doing.


%% Listen to the filtered sound
% sound(y_down_sampled, SampleRate_down_sampled);

%% Save the filtered signal to a wave file
wavwrite(y_down_sampled/(0.1+max(y_down_sampled)), SampleRate_down_sampled, NbitsPerSample, 'voice_record_down_sampled.wav');





