function [reverbsignal] = SpringReverbSimulator(signal, fs, impulseresponse, npoles)
%
% SpringReverbSimulator.m 
%
% Analyses spring reverb system impulse responses to construct a filter-based 
% representation of the specified system in conjunction with a fixed algorithmic 
% reverberator. The function also processes input signal for immediate playback, 
% wave file writing and comparative plotting.
%
%
% Syntax:
% [reverbsignal] = ReverbSimulator(signal, fs, impulseresponse, npoles)
%
%       signal              =   Signal to be processed
%       fs                  =   Sample rate in Hertz
%       impulseresponse     =   Impulse response of system to be simulated
%       npoles              =   Number of pole coefficients for LPC stage
%       reverbsignal        =   Reverberant signal output
%
%              
% Construction:                 
%                                           FILTER 1
%                                           FILTER 2
% INPUT -> APC ANALYSIS -> APC FILTER ->    FILTER 3    ->  AP1  ->  AP2  ->  AP3  -> AP4   -> OUTPUT
%                                           FILTER 4
%                                           FILETR 5
%
% Process:
% From the initial colouration synthesised by the LPC-based filter stage,
% the reverberation is augmented by a fixed reverberator constructed from 
% both parallel and series delay-line filters. The fixed reverberator is
% based on generalised and typical spring reverb properties.
%
% Example:
% Simulation of Roland RE-201 Space Echo spring reverb unit. Preliminary
% steps are as follows.
%
% [impulseresponse, fs, nbits] = wavread(RolandRE201ImpulseResponse.wav);
% signal = wavread(WaveFileToProcess.wav);
% npoles = 12;
% reverbsignal = SpringReverbSimulator(signal, fs, impulseresponse, npoles);
%


% LPC ANALYSIS

polecoeff = lpc(impulseresponse, npoles);
impresps = impz(1, polecoeff, 256, fs); 
impresps = impresps./max(abs(impresps));
hlarge = freqz(1, polecoeff, 1024, fs);
lpcmag = mag2db(abs(hlarge));
magoffset = median(lpcmag);
lpcmag = lpcmag-magoffset;
       
% LPC Filter Coefficients Applied

treatedsignal = filter(impresps, lpcmag, signal); 
treatedsignal = treatedsignal./max(abs(treatedsignal));



% PARALLEL FILTER BANK

% Parallel Filter 1

d = 2106;
gain = 0.7;
   
    b = [gain zeros(1, d) -1];
    a = [1 zeros(1, d) -gain];
    pfilt1 = filter(b, a, treatedsignal);


% Parallel Filter 2

d = 3804;
gain = 0.7;
   
    b=[gain zeros(1,round(d)) -1];
    a=[1 zeros(1,round(d)) -gain];
    pfilt2 = filter(b, a, treatedsignal);


% Parallel Filter 3

d = 5440;
gain = 0.7;
   
    b=[gain zeros(1,round(d)) -1];
    a=[1 zeros(1,round(d)) -gain];
    pfilt3 = filter(b, a, treatedsignal);


% Parallel Filter 4

d = 7880;
gain = 0.7;
   
    b=[gain zeros(1,round(d)) -1];
    a=[1 zeros(1,round(d)) -gain];
    pfilt4 = filter(b, a, treatedsignal);


% Parallel Filter 5

d = 8040;
gain = 0.7;
   
    b=[gain zeros(1,round(d)) -1];
    a=[1 zeros(1,round(d)) -gain];
    pfilt5 = filter(b, a, treatedsignal);


% Sum of Parallel Filters

sumsignal = pfilt1 + pfilt2 + pfilt3 + pfilt4 + pfilt5;



% SERIES FILTERS

% All-pass Filter 1

delay = 220;
gain = 0.7;

b = [-gain zeros(1, delay-1) 1];
a = [1 zeros(1, delay-1) -gain];
ap1signal = filter(b, a, sumsignal);

% All-pass Filter 2

delay = 441;
gain = 0.5;

b = [-gain zeros(1, delay-1) 1];
a = [1 zeros(1, delay-1) -gain];
ap2signal = filter(b, a, ap1signal);

% All-pass Filter 3

delay = 713;
gain = 0.3;

b = [-gain zeros(1, delay-1) 1];
a = [1 zeros(1, delay-1) -gain];
ap3signal = filter(b, a, ap2signal);

% All-pass Filter 4

delay = 1124;
gain = 0.1;

b = [-gain zeros(1, delay-1) 1];
a = [1 zeros(1, delay-1) -gain];
reverbsignal = filter(b, a, ap3signal);



% OUTPUT STAGE

% Normalise Output

reverbsignal = reverbsignal./max(abs(reverbsignal));

% Write WAVE File

wavwrite(reverbsignal, fs, 'SpringReverbSimulatorOUTPUT.wav');

% Playback Reverberant Signal

sound(reverbsignal, fs);

% Plot Reverberant Signal and Original Signal

subplot(2,1,1), plot(reverbsignal)
title ('Reverberant Signal')
subplot(2,1,2), plot(signal)
title ('Original Signal')

end

