// Signals for psychoacoustic experiments

caso = 6;

select caso
case 0
   // 0. Bessel functions
   m = (0:0.1:10)'
   x = besselj([0,1,2,3,4],m);
   clf
   plot(m,x)
   title("Bessel functions")
   legend("J0", "J1", "J2", "J3", "J4");

case 1
   // Sine wave
   // Duration in s
   T = 2;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';   
   // Frequency in Hz
   f = 1000; 
   // Signal
   x = sin(2*%pi*f*t);
   // Play sound
   playsnd(0.1*x, Fs, 1)
   // Signal spectrum
   N = 4096
   F = (0:N/2-1)'*Fs/N;
   X = fft(x(1:N),-1);
   X = 2/N * abs(X(1:N/2));
   // Plot of signal and its spectrum 
   clf
   subplot(2,1,1)
   plot(t(1:220),x(1:220))
   subplot(2,1,2)
   plot("ll",F(1:200),X(1:200))

case 2
   // Polyharmonic
   // Duration in s
   T = 5;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';   
   // Number of harmonics
   n = 5;    
   // Frequencies in Hz
   f = 1000 * (1:n);         
   // Amplitudes
   A = 0.5 * 1 ./ (1:n);   
   // Signal initialization
   x = zeros(t);            
   for k=1:n
      // Signal update
      x = x + A(k)*sin(2*%pi*f(k)*t);
   end
   // Play sound
   playsnd(0.1*x(1:2*Fs), Fs, 1)
   // Signal spectrum
   N = 4096
   F = (0:N/2-1)'*Fs/N;
   X = fft(x(1:N),-1);
   X = 2/N * abs(X(1:N/2));
   clf
   subplot(2,1,1)
   plot(t,x) 
   subplot(2,1,2)
   plot(F,X)

case 3
   // Inharmonic
   T = 2;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';   
   // Frequencies in Hz
   f = [1000, 1245, 1722.8, 2325] 
   // Number of partials
   n = length(f);        
   // Amplitudes
   A = [0.15 0.3 0.1 0.05];   
   // Signal initialization
   x = zeros(t);            
   for k=1:n
      // Signal update
      x = x + A(k)*sin(2*%pi*f(k)*t);
   end
   // Application of a decreasing exponential envelope
   // to simulate a bell
   xx = exp(-t/0.7).*x;
   // Play sounds (constant and decreaing amplitude)
   playsnd(0.2*x, Fs, 1)
   playsnd(0.2*xx, Fs, 1)
   // Signal spectrum
   // FFT size
   N = 4096;
   // Frequency vector
   F = (0:N/2-1)'*Fs/N;
   // Signal FFT
   X = fft(x(1:N), -1);
   // Conversion to real spectrum
   X = 2/N * abs(X(1:N/2));
   // Periodic Hann window
   Hann = window('hn', N+1)';
   Hann = Hann(1:N);
   // Energy attenuation due to Hann window
   Q = sqrt(mean(Hann.^2))
   // Hann-windowed spectrum
   X1 = fft(x(1:N).*Hann,-1); 
   X1 = 2/N * abs(X1(1:N/2));

   clf
   subplot(2,1,1)
   plot(t,x) 
   title("Inharmonic signal") 
   xlabel("t  (s)")
   ylabel("x")
   gca().data_bounds = [0, -0.7; 0.04, 0.7];  
   subplot(2,1,2)
   plot(F,[X,X1/Q])
   xlabel("f  (Hz)")
   ylabel("|X|")
   title("Spectrum of inharmonic signal") 
   legend("No window", "Hann window");
   gca().data_bounds = [0, 0; 5000, 1.1*max(X)];  

case 4
   // White noise with gaussian distribution
   // Duration in s
   T = 20;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';
   x = grand(length(t),1,"nor",0,0.25);
   N = 4096;
   // Play sound
   playsnd(0.1*x(1:2*Fs), Fs, 1)
   // One-frame fft spectrum
   F = (0:N/2-1)'*Fs/N;
   X = fft(x(1:N));
   X = 2/N * abs(X);
   // Average fft spectrum 
   M = floor(length(x)/N);
   x1 = x(1:M*N);
   x1 = matrix(x1,N,M);
   // Column-wise FFT
   X1 = fft(x1,-1,1);
   // Energy average of columns and amplitude adjustmnt
   X1 = 2/N * sqrt(mean(abs(X1).^2,"c"));    
   clf
   subplot(2,1,1)
   plot(t,x)
   subplot(2,1,2)
   plot(F,X1(1:N/2))
   plot(0,0)

case 5
   // White noise with uniform distribution
   // Duration in s
   T = 2;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';
   x = grand(length(t),1,"unf",-1,1);
   N = 8192;
   F = (0:N/2-1)'*Fs/N;
   // Play sound
   playsnd(0.1*x, Fs, 1)    
   // Single window FFT spectrum
   X = fft(x(1:N));
   X = 2/N * abs(X);

   clf
   subplot(2,1,1)
   plot(t,x)
   subplot(2,1,2)
   plot(F,X(1:N/2)) 

case 6
   // Pink noise
   // Duration in s
   T = 100;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';
   // Pink noise
   // N for pink noise generation
   N = 16384;    
   exec pinknoise.sci;
   x = pinknoise(T, Fs, N, 1);
   // Play sound
   playsnd(0.1*x(1:5*Fs), Fs, 1)    
   // N for spectrum analysis
   N = 32768;
   // Single-window FFT spectrum
   // Frequency vector
   F = (0:N/2-1)'*Fs/N;       
   // FFT of first N samples
   X = fft(x(1:N));            
   // Adjust FFT for physical amplitudes
   X = 2/N * abs(X(1:N/2));   
   // Average FFT spectrum 
   // Number of complete windows
   M = floor(length(x)/N);    
   // Limit x to complete windows
   x1 = x(1:M*N);             
   // Convert x1 to matrix, one column per window
   x1 = matrix(x1, N, M);     
   // Column-wise FFT
   X1 = fft(x1, -1, 1);
   // Adjust FFT for physical amplitudes
   X1 = 2/N * abs(X1(1:N/2,:));   
   // Energy average of columns 
   X1 = sqrt(mean(X1.^2,"c"));
   // Linear regression to fit FFT spectrum except for the 50 lowest 
   // frequencies to avoid the effect of low-frecuency anomalies  
   [a, b] = reglin(log(F(50:$)'),log(X1(50:$)'))
   X2 = exp(a*log(F) + b);
   // Signal and spectrum plot    
   clf
   subplot(2,1,1)
   plot2d('nn', t(1:Fs),x(1:Fs))
   // Axes for subplot 1
   sub1 = gca();
   subplot(2,1,2)
   plot2d('ll', F, [X1,X2])
   // Axes for subplot 2
   sub2 = gca();
   // Select colors
   sub1.children.children.foreground = 1;
   sub2.children.children(1).foreground = 1;
   sub2.children.children(2).foreground = 2;

case 7
   // Brownian noise
   // Duration in s
   T = 2;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';
   // Create brownian noise
   exec browniannoise.sci;
   x = browniannoise(T, Fs);
   // Brownian noise playback
   playsnd(0.1*x,Fs)
   // N for spectrum analysis
   N = 8192;
   F = (0:N/2-1)'*Fs/N;
   X = fft(x(1:N));
   X = 2/N * abs(X(1:N/2));
   clf
   subplot(2,1,1)
   plot(t,x)
   subplot(2,1,2)
   plot("ll", F, X) 
   gca().grid = [1 1]

case 8
   // Band noise
   // Duration in s
   T = 5;
   // Sample rate
   Fa = 44100;
   // Center frequency
   fo = 1000;
   f1 = fo*2^(-1/6);
   f2 = fo*2^(1/6);
   // Size of the band noise filtering window 
   N = 8192;
   // Effective value
   Xef = 1
   // Create band noise
   exec bandnoise.sci;
   x = bandnoise(f1, f2, T, Fs, Xef, N);
   // Time vector
   t = (0:length(x)-1)/Fs;
   // Play band noise
   playsnd(x, Fs, 1)
   // N for spectrum analysis
   N = 8192;
   // Frequency vector
   F = (0:N/2-1)'*Fs/N;
   // FFT spectrum
   X = fft(x(1:N));
   X = 2/N * abs(X(1:N/2));
   // Mean FFT spectrum
   exec fftmean.sci;
   w = "hann";
   L = N;
   [y2, y, ymin, ymax, ff] = fftmean(x, Fs, N, w, L)
   
   // Plot signal and spectrum
   clf();
   subplot(3,1,1)
   plot(t,x)
   xlabel("t  (s)")
   ylabel("x(t)")    
   subplot(3,1,2)
   plot("ln", F, X)
   xlabel("f  (Hz)")
   ylabel("|X(f)|") 
   gca().grid = [1 1];
   gca().data_bounds = [fo/2 0; fo*2, Xef/2];
   subplot(3,1,3)
   plot("ln", F, y2)
   xlabel("f  (Hz)")
   ylabel("|X(f)|") 
   gca().grid = [1 1];
   gca().data_bounds = [fo/2 0; fo*2, Xef/2];
   
case 9
   // Logon
   // Duration in s
   T = 2;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';      // Frequency in Hz
   f = 2000;
   // Sigma in s 
   sigma = 0.001;
   // Gaussian envelope 
   env = 1/sqrt(2*%pi)/sigma * exp(-((t-T/2)/sigma).^2);
   x = env .* sin(2*%pi*f*t);
   N = 8192;
   F = (0:N/2-1)'*Fs/N;
   TN = T*Fs/2 + (-N/2:N/2-1)';
   X = fft(x(TN));
   X = 2/N * abs(X);
   clf
   subplot(2,1,1)
   TT = floor(T*Fs/2-1000:T*Fs/2+1000)';
   plot(t(TT),x(TT),t(TT),env(TT))    
   subplot(2,1,2)
   plot(F(1:1000),X(1:1000)) 
   playsnd(x,Fs)
end


N = 16
q = window("hn",N+1);
q = q(1:N);
q1 = [q, q, q, q];
q2 = [q(N/2+1:N), q, q, q, q(1:N/2)];
Q = q1 + q2