// Ejemplo de señal con aliasing
//
// Para visualizar el problema del aliasing se generan
// cuatro ondas: Senoidal, cuadrada sin aliasing, cuadrada 
// con poco aliasing y cuadrada con mucho aliasing. La 
// cuadrada sin aliasing se genera sumando armónicos hasta
// Fs/2, La cuadrada con poco aliasing, sumando armónicos
// que exceden la tasa de Nyquist Fs/2. La cuadrada con 
// mucho aliasing tomando el signo de una senoide.
// 
// Para evitar contaminar el ejemplo se usó una Señal sin
// derrame espectral. Para ello se eligió la frecuencia
// de modo que quepa un número entero de periodos en la
// ventana de análisis:
//
// fo = m*Fs/N
//
// -----------------------------------------------------
// Federico Miyara
// 2024-10-24

clear
// Tasa de muestreo
Fs = 44100
// Tamaño de la ventana para FFT
N = 4096

// Número de línea espectral donde se encontrará
// la frecuencia fundamental (armónico 1)
m = 95
// Periodo de las ondas para que no haya derrame
// espectral (se consigue haciendo caber exactamente
// m periodos.
// NOTA: salvo que además m sea una potencia de 2,
// la señal digital obtenida no es periódica, aunque
// sí cumple el teorema de muestreo (no hay alias) 
// representa una onda periódica en tiempo continuo 
To = N/Fs/m
// Frecuencia fundamental
fo = 1/To;
// Duración de las señales 
T = N/Fs;
// Vector de tiempos
t = (0:N-1)/Fs;
// Señal senoidal
x = sin(2*%pi*fo*t);
// Onda cuadrada sin alias
y = 0;
for k=1:2:floor(Fs/2/fo)
   y = y + 4/%pi * sin(2*%pi*k*fo*t) / k;
end

// Onda cuadrada con alias a partir de su 
// serie de Fourier con 35 armónicos
z = 0;
for k=1:2:35
   z = z + 4/%pi * sin(2*%pi*k*fo*t) / k;
end

// Onda cuadrada con alias generada 
// aplicando la función signo a una
// senoide ligeramente desfasada para
// evitar valores 0. 
w = sign(sin(2*%pi*fo*t + 0.0001));

// Espectros
X = abs(fft(x(1:N)))*2/N;
Y = abs(fft(y(1:N)))*2/N;
Z = abs(fft(z(1:N)))*2/N;
W = abs(fft(w(1:N)))*2/N;
W = max(0.000001, W);

// Vector de frecuencias
f = (0:N/2-1)*Fs/N;

caso = "log"
select caso
case "lin"
   scf(1);
   clf(1);
   subplot(4,2,1)
   plot(t(1:2*N/m), x(1:2*N/m));
   gca().data_bounds = [0 -1.2; t(2*N/m) 1.2];
   subplot(4,2,2)
   plot(f, X(1:N/2))
   gca().data_bounds = [0 0; Fs/2 4/%pi];
   subplot(4,2,3)
   plot(t(1:2*N/m), y(1:2*N/m));
   gca().data_bounds = [0 -1.2; t(2*N/m) 1.2];
   subplot(4,2,4)
   plot(f, Y(1:N/2))
   gca().data_bounds = [0 0; Fs/2 4/%pi];
   subplot(4,2,5)
   plot(t(1:2*N/m), z(1:2*N/m));
   gca().data_bounds = [0 -1.2; t(2*N/m) 1.2];
   subplot(4,2,6)
   plot(f, Z(1:N/2))
   gca().data_bounds = [0 0; Fs/2 4/%pi];
   subplot(4,2,7)
   plot(t(1:2*N/m), w(1:2*N/m));
   gca().data_bounds = [0 -1.2; t(2*N/m) 1.2];
   subplot(4,2,8)
   plot(f, W(1:N/2))
   gca().data_bounds = [0 0; Fs/2 4/%pi];
case "log"
   scf(1);
   clf(1);
   subplot(4,2,1)
   plot(t(1:2*N/m), x(1:2*N/m));
   gca().data_bounds = [0 -1.2; t(2*N/m) 1.2];
   subplot(4,2,2)
   plot("nl", f, X(1:N/2))
   gca().data_bounds = [0 0.001; Fs/2 4/%pi];
   subplot(4,2,3)
   plot(t(1:2*N/m), y(1:2*N/m));
   gca().data_bounds = [0 -1.2; t(2*N/m) 1.2];
   subplot(4,2,4)
   plot("nl", f, Y(1:N/2))
   gca().data_bounds = [0 0.001; Fs/2 4/%pi];
   subplot(4,2,5)
   plot(t(1:2*N/m), z(1:2*N/m));
   gca().data_bounds = [0 -1.2; t(2*N/m) 1.2];
   subplot(4,2,6)
   plot("nl", f, Z(1:N/2))
   gca().data_bounds = [0 0.001; Fs/2 4/%pi];
   subplot(4,2,7)
   plot(t(1:2*N/m), w(1:2*N/m));
   gca().data_bounds = [0 -1.2; t(2*N/m) 1.2];
   subplot(4,2,8)
   plot("nl", f, W(1:N/2))
   gca().data_bounds = [0 0.001; Fs/2 4/%pi];
end

