// simulation of imupulse response measurement using continuous TSP (Time Stretched Signal)
// http://www.cepstrum.co.jp/

clear;                                        // clear all variable
stacksize('max');                             // allocate maximum working memory
rand('uniform');                              // use uniform random noise

len=16384;                                    // TSP signal period
alpha=-1.0;                                   // coefficient to generate continuous TSP
noise_amp=0.1;                                // noise level

// generate continuous TSP signal in frequency domain
x=zeros(1:len);
for i=1:len/2
  c=2*%pi/len*alpha*i^2;
  x(i+1)=cos(c)+%i*sin(c);
end;
for i=1:len/2-1
  x(len+1-i)=conj(x(i+1));
end;

// generate time domain continuous TSP signal
tsp=real(fft(x, 1));                          // IFFT

revtsp=mtlb_fliplr(tsp);                      // generate reverse TSP signal

gain_factor=sqrt(2*len)/2;                    // calculate scaling factor
tsp=tsp*gain_factor;                          // scaling 

tsp2=[tsp, tsp];                              // generate 2 period of TSP signal

// generate impulse response example for simulation
h=(2*rand(1:len)-1).*(0.999^(1:len));
h=[zeros(1:100), h(1:len-100)];

res=convol(h, tsp2);                          // calculate response signal
res=res+(2*rand(1:length(res))-1)*noise_amp;  // add noise  
res=res(len+1:len+len);                       // extract response signal

// calculate frequency response from response signal
ampspec1=20*log10(abs(fft(res/gain_factor, -1))+1e-19);  

// calculate impulse response 
imp=fft(fft(res, -1).*fft(revtsp, -1), 1);    // circular convolution by FFT/IFFT
imp=imp/gain_factor;                          // scaling 

// calculate energy decay curve by Schroeder's method
decay=zeros(1:len);
engsum=0;
for i=len:-1:1
  engsum=engsum+imp(i)*imp(i);
  decay(i)=engsum;
end
decay=decay/decay(1);
decay=10*log10(decay);                        // calculate log curve of energy decay

// calculate frequency response from measured impulse response
ampspec2=20*log10(abs(fft(imp, -1))+1e-19);  

scf(0); clf; plot(tsp2);                      // plot 2 period of continuous TSP signal
scf(1); clf; plot(h);                         // plot impulse response to be measured
scf(2); clf; plot(res);                       // plot response signal of the system
scf(3); clf; plot(imp);                       // plot measured impulse response
scf(4); clf; plot(decay);                     // plot energy decay

// plot frquency response 
scf(5); clf; plot2d(ampspec1, rect=[0, 0, len/2, 40], style=2);       
scf(6); clf; plot2d(ampspec2, rect=[0, 0, len/2, 40], style=5);    

savewave('tsp.wav', 0.9*tsp2, 8000);          // save TSP continuous TSP singal (Fs=8[kHz])