%% SystemResponseToolbox Example Script
clear all; close all; clc;

%% INPUT PARAMETERS
fs = 48000;
N_fft = 4096;
InputSignal = wavread('whitenoise_3s_48kHz.wav');
OutputSignal = wavread('mic_134.wav');% This is taken from the Soundbender measurement session of the 25102012
OutputSignal = OutputSignal(2*fs:5*fs);% Selecting the signal of the first loudspeaker of the Soundbender
EstimatorType = 'H1';

%% PROCESSING STAGE

% CHOICE 1: get the results of the estimation as single variables (that one
% can use in the "traditional way" AKA "functions and variables way"...)
[h, H, COH, f, t] = SystemResponseToolbox.HEstimator(InputSignal, OutputSignal, N_fft, fs, EstimatorType); % Estimating the IR and FRF via H1 Estimator
% The above command is enough if you want to use the toolboox and integrate
% it with your code. In this case, the output of the HEstimator is enough.

% Let's now have a look to the Object Oriented methods.

% CHOICE 2: get the results in an Object (instead of variables): "The new way"!
SysRespObj = SystemResponseToolbox(InputSignal, OutputSignal, N_fft, fs, EstimatorType, 'GetTheMinimumPhase');
% In few words, all the information contained in the h, H, COH, f, and t
% variables that we normally use in our script are contained in an Object,
% named SysRespObj (or you can give to it the name you prefer) and they can
% be accessed with this command
% >> SysRespObj.IRcoefficients %this gives the same values as the variable
% h of CHOICE1

% List the Properties/Methods of the Object, just to get an idea (you do not always have to do it)...
SysRespObj
methods(SysRespObj) %The methods are "functions" that can help us in performing some
%operation that are often repeated for a given dataset (i.e., plotting,
%etc.)

%% Plotting

% By the following, a list of examples of the possibilities offered by the
% Object Paradigm for  plotting

% Plots the results of the estimation (IR, FRF, COH) in the same plot. This
% could be useful to check the results of our estimation
SystemResponseToolbox.Plot_Estimation_Results(SysRespObj);

% Or you may want to plot the single bits of the System Response, for
% example the Impulse response of our System
figure; SysRespObj.PlotIR(SysRespObj);

% Or you may want to combine the figures in various subplots, such as:
figure; 
subplot(1,2,1); SysRespObj.Plot_Magnitude_FRF(SysRespObj); grid on;
subplot(1,2,2); SysRespObj.Plot_Phase_FRF(SysRespObj); grid on;

% If you have two Objects (named, for example, SysRespObj1 and SysRespObj2)
% you can plot them on the same figure, i.e.
% >> figure; 
% >> SysRespObj1.Plot_Magnitude_FRF(SysRespObj1);
% >> hold all; SysRespObj2.Plot_Phase_FRF(SysRespObj2);

break;
% Log-smoothing of the estimated IR
%SystemRespObj.IRcoefficients = sequence_log_smoothing(SystemRespObj.IRcoefficients, SystemRespObj.fs, 1);


% %% Minimum and Maximum Phase
% 
% % Minimum phase
% h_minimum_phase = minimumphasesequence(h);
% IRminimphase = ImpulseResponse(h_minimum_phase);
% H_min_phase = fft(h_minimum_phase); H_min_phase = H_min_phase(1:N_fft/2 + 1);
% FRFminimphase = TransferFunction(H_min_phase, f, fs, N_fft);% Creates a new TransferFunction object with vector H
% 
% % Maximum Phase
% H_max_phase = H./H_min_phase;
% h_max_phase = ifft(H_max_phase, N_fft, 'symmetric');% Creates a new IR object with coeffients h
% IRmaxphase = ImpulseResponse(h_max_phase);
% FRFmaxphase = TransferFunction(H_max_phase, f, fs, N_fft);% Creates a new TransferFunction object with vector H
% 
% %% Inverse Filtering
% 
% % Calculating the ORIGINAL Inverse Filter
% InverseFilterFRF = 1./H;
% InverseFilter = TransferFunction(InverseFilterFRF, f, fs, N_fft);% Creates a new TransferFunction object with vector H
% 
% % Inverse filter of the minumum phase component
% InverseFilterMinPhaseFRF = 1./H_min_phase;
% %InverseFilterMinPhaseFRF = InverseFilterMinPhaseFRF.*exp(-1i*(pi - deg2rad(0))*transpose(f));
% InverseFilterMinPhase = TransferFunction(InverseFilterMinPhaseFRF, f, fs, N_fft);% Creates a new TransferFunction object with vector H

%% Plots

% Plot the INPUT/OUTPUT signals
% figure; subplot(1,2,1);
% plot(InputSignal); title('Input Signal');
% subplot(1,2,2); plot(OutputSignal); title('Output Signal');

%IR
figure; IR.PlotIR(IR);
hold all; IRminimphase.PlotIR(IRminimphase);
hold all; IRmaxphase.PlotIR(IRmaxphase);
convolution = conv(IRminimphase.coefficients,IRmaxphase.coefficients);
IRConvolution = ImpulseResponse(convolution);
hold all; IR.PlotIR(IRConvolution);
legend('Original IR', 'IR (Min Phase component)', 'IR (Max Phase component)', 'Convolution between Min and Max phase');

% FRF
figure; 
subplot(1,2,1); FRFminimphase.Plot_Magnitude_FRF(FRFminimphase);
hold all; FRF.Plot_Magnitude_FRF(FRF); grid on;
subplot(1,2,2); FRFminimphase.Plot_Phase_FRF(FRFminimphase, 1);  
hold all; FRF.Plot_Phase_FRF(FRF); grid on;

figure;
subplot(1,2,1); InverseFilter.Plot_Magnitude_FRF(InverseFilterMinPhase);
hold all; InverseFilterMinPhase.Plot_Magnitude_FRF(InverseFilter); grid on;
title('Magnitude');

subplot(1,2,2); InverseFilter.Plot_Phase_FRF(InverseFilterMinPhase);
hold all; InverseFilterMinPhase.Plot_Phase_FRF(InverseFilter); grid on;
title('Phase')
legend('Inverse filter Minimum phase component','Original inverse filter');