%% load wave file 
% wave files have inside information about
% sampling rate, usually left and right channel amplitudes, how many bits per sample etc. 
[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 often 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);
%___%
%___%%% lets find out the sample rate
%___%dt_per_point=diff(t);  % notice that matlab has this shortcut dt(i)=t(i+1)-t(i)
%___%
%___%% let's 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
%___%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
% 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 frequency at the right we need to find index of sorted freq
[dummy,index]=sort(spectrum_freq);
plot(spectrum_freq(index), abs(Y(index)) ); % notice abs generaly Y is complex
xlabel('Frequency (Hz)');
ylabel('FFT amplitude');
title('Spectrum amplitudes');
% have a look at this plot and chose proper filter roll off frequency

%% let's low pass data
Freq3db=100;
[spectrum_filtered,g_vs_freq]=freq_space_low_pass_filter(Y, SampleRate, Freq3db);

%% filter frequency gain 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');

%% filtered and unfiltered spectrum
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');

%% construct filtered signal
y_filtered=ifft(spectrum_filtered);

%% check that filter worked right so 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 was good
y_filtered=real(y_filtered);

%% compare filtered and raw signals
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');

%% output filtered sound
% sound(y_filtered, SampleRate);

%% save the wave file with filtered sound
wavwrite(y_filtered, SampleRate, NbitsPerSample, 'voice_record_filtered.wav');


%% down sample
% since a lot of amplitudes in our spectrum are zero
% we will make truncated spectrum with frequencies only with in filter band
% the band is depends on filter for RC type 10*Freq3db is a good estimate
% Remember!  Down sampling a not filtered spectrum is a waste of data points.
% You can do it by simple decimation (just skip take only Mth point)!
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;
plot(t, y_filtered,'r-');
plot(t_down_sampled, y_down_sampled,'m-');
title('Signal vs time');
xlabel('Time (S)');
ylabel('Amplitude');
legend('raw', 'filtered', 'down sampled');
% you probably will not able to distinguish filtered and down sampled version


%% output filtered sound
% sound(y_down_sampled, SampleRate_down_sampled);

%% save the wave file with filtered sound
wavwrite(y_down_sampled, SampleRate_down_sampled, NbitsPerSample, 'voice_record_down_sampled.wav');





