
// Signals for psychoacoustic experiments
//
// Some examples (not quite organized) of functions 
// and signals frequently used in psychoacoustics
//
// Set caso = 1, ..., 9 to see and listsn to different
// signals.
//
// 0: Not really a signal, but Bessel functions which
//    are related to FM spectrrum
// 1: Sine waves
// 2: Polyharmonic signals
// 3: Inhaermonic signals 
// 4: Gaussian white noise
// 5: Uniform whiye noise
// 6: Pink noise
// 7: Brownian noise
// 8: Band noise
// 9: Logon (gaussian-envelope beep
//
// ---------------------------
// Federico Miyara
// 2019-05-14
// 2024-10-24

clear
close(winsid());

caso = 6;

select caso
case 0
   // 0. Bessel functions
   m = (0:0.1:10)'
   x = besselj([0,1,2,3,4],m);
   scf(0)
   clf(0)
   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));
   // 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))/Q;
   
   // Plot of signal and its spectrum 
   scf(1);
   clf(1);
   subplot(2,1,1)
   plot(t(1:220),x(1:220))
   subplot(2,1,2)
   plot("ll", F(1:200), X(1:200), F(1:200), X1(1:200))
   legend("No window", "Hann window",2)

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));
   // 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))/Q;
   
   scf(2);
   clf(2);
   subplot(2,1,1)
   plot(t,x) 
   gca().data_bounds = [0, -0.7; 0.005, 0.7];  
   subplot(2,1,2)
   plot(F, X, F, X1)
   legend("No window", "Hann window")

case 3
   // Inharmonic
   T = 4;
   // 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(1:2*Fs), 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));

   scf(3);
   clf(3)
   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 2 s of 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 
   // Number of frames
   M = floor(length(x)/N);
   // Limit signal to an integer number of frames
   x1 = x(1:M*N);
   // Resize as M columns
   x1 = matrix(x1,N,M);
   // Column-wise FFT
   X1 = 2/N * abs(fft(x1,-1,1));
   // Energy average of columns and amplitude adjustmnt
   X1avg = sqrt(mean(X1.^2,"c"));    
   
   scf(4);
   clf(4);
   subplot(3,1,1)
   plot(t,x)
   xlabel("t  (s)")
   ylabel("x(t)")
   subplot(3,1,2)
   plot(F, X(1:N/2))
   xlabel("f  (Hz)")
   ylabel("|X|_average")
   gca().data_bounds = [0 0; 20000 max(X)];
   subplot(3,1,3)
   plot(F, X1avg(1:N/2))
   xlabel("f  (Hz)")
   ylabel("|X|_average")
   gca().data_bounds = [0 0; 20000 max(X)];

case 5
   // White noise with uniform distribution
   // Duration in s
   T = 20;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector
   t = (0:1/Fs:T)';
   x = grand(length(t),1,"unf",-1,1);
   N = 4096;
   F = (0:N/2-1)'*Fs/N;
   // Play sound
   playsnd(0.1*x(1:2*Fs), Fs, 1)    
   // Single window FFT spectrum
   X = fft(x(1:N));
   X = 2/N * abs(X);
   // Average fft spectrum 
   // Number of frames
   M = floor(length(x)/N);
   // Limit signal to an integer number of frames
   x1 = x(1:M*N);
   // Resize as M columns
   x1 = matrix(x1,N,M);
   // Column-wise FFT
   X1 = 2/N * abs(fft(x1,-1,1));
   // Energy average of columns and amplitude adjustmnt
   X1avg = sqrt(mean(X1.^2,"c"));    
   
   scf(5);
   clf(5);
   subplot(3,1,1)
   plot(t, x)
   xlabel("t  (s)")
   ylabel("x(t)")   
   gca().data_bounds = [0, -1.2; 1, 1.2];
   subplot(3,1,2)
   plot(F, X(1:N/2)) 
   xlabel("f  (Hz)")
   ylabel("|X|")
   gca().data_bounds = [0 0; 20000 max(X)];
   subplot(3,1,3)
   plot(F, X1avg(1:N/2))
   xlabel("f  (Hz)")
   ylabel("|X|_average")
   gca().data_bounds = [0 0; 20000 max(X)];
   
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    
   scf(6);
   clf(6);
   subplot(3,1,1)
   plot(t(1:Fs), x(1:Fs))
   xlabel("t  (s)")
   ylabel("x(t)")
   subplot(3,1,2)
   plot("ll", F, X)
   xlabel("f  (Hz)")
   ylabel("|X|")
   subplot(3,1,3)
   plot("ll", F, X1)
   plot("ll", F, X2, "r")
   xlabel("f  (Hz)")
   ylabel("|X|_med, |X|_teor") 
   legend("|X|_med", "|X|_teor", 3)  

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));
   scf(7);
   clf(7);
   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 = 50;
   // Sample rate
   Fs = 44100;
   // Center frequency
   fo = 1000;
   // Octave fraction
   b = 1/12
   f1 = fo*2^(-b/2);
   f2 = fo*2^(b/2);
   // 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(0.2*x(1:3*Fs), 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
   scf(8);
   clf(8);
   subplot(3,1,1)
   plot(t,x)
   xlabel("t  (s)")
   ylabel("x(t)")    
   gca().data_bounds = [0 -max(abs(x))*1.1; 0.5, max(abs(x))*1.1];
   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, max(X)];
   subplot(3,1,3)
   plot("ln", F, y2)
   xlabel("f  (Hz)")
   ylabel("|X(f)|_med") 
   gca().grid = [1 1];
   gca().data_bounds = [fo/2 0; fo*2, max(X)];
   
case 9
   // Logon
   // Duration in s
   T = 2;
   // Sample rate in Hz
   Fs = 44100;
   // Time vector in s
   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);
   
   // Play sound
   playsnd(x,Fs)

   scf(9);
   clf(9);
   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)) 
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



