
// Testing the function spectrogram()

clear

// Sample rate i Hz
Fs = 44100
// Size of FFT window 
N = 1024
// Ofset in samples between FFT windows
M = 256
// Tapering window function
w = "hann"
// Type of graphic output
grph = "logcmap";

// Time vector in s
t = 0:1/Fs:2;
// Carrier frequency in Hz
fo = 8000
// Modulation frequency in Hz
fm = 0.5
// Frequency deviation in Hz
df = 4000
// Frequency vector
f = fo + df*sin(2*%pi*fm*t);
// Phase computed by numerical integration 
phi = 2*%pi*cumsum(f)/Fs;

// Signal generation
caso = 5
select caso
case 1
    // Two harmonics  
    x = 0.4*sin(phi/2) + 0.3*sin(phi) + 0.1*sin(3/2*phi);
    grph = "logcmap";
case 2
    // Three harmonics with aliasing 
    x = 0.4*sin(phi) + 0.3*sin(2*phi);
case 3
    // Two harmonics with aliasing and distortion
    // by truncation
    x = sin(phi) + sin(2*phi);  
    x = min(x, ones(x));  
case 4
    Fs = 22050;
    t = 0:1/Fs:1*(1-%eps);
    x = 0.4 + [1.5*sin(2*%pi*800*t) sin(2*%pi*1200*t)]; 
    w = "flattop"
    grph = "lincmap";
 case 5
    ruta = "E:\work_scilab\";
    archivo = "tono corto para espectrograma.wav";
    w = "hann"
    [x, Fs] = wavread(ruta+ archivo);
    grph = "logcmap"         
end

// Save wave as a wav file to compare with external
// software such as Audacity
//wavwrite(x,"d:\hola\fmod.wav")

exec spectrogram.sci;

close()
z = spectrogram(x, N, M, Fs, w, grph, "jet");
xlabel("t  [s]")
ylabel("f  [Hz]")



