[Scilab-users] spectrogram

[Federico Miyara]

Dear All,

I noticed there is no spectrogram function in Scilab that provides both 
graphic and numerical output. There is a function mapsound(), which 
plots a spectrogram without numerical information and without control 
over important aspects such has FFT length or windowing function. By the 
way, in the field of digital sound processing the usual name for this is 
spectrogram, not mapsound, which is a bit confusing.

I created a spectrogram() function. Even if it doesn't yet comply with 
the recommended style for Scilab functions (particularly regarding error 
handling or manual page), I think it might be useful for some users. I 
attach the function and a test script.

Regards,

Federico Miyara




-- 
El software de antivirus Avast ha analizado este correo electrónico en busca de virus.
https://www.avast.com/antivirus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20210711/f91d090b/attachment.htm>
-------------- next part --------------
function [z, t, f] = spectrogram(x, N, Fs, M, w, grph)
    // Spectrogram of a signal
    //
    // Usage
    //       [z, t, f] = spectrogram(x, N, Fs, M, w, grph)
    // where
    //       x:    signal vector
    //       Fs:   sample rate in Hz
    //       N:    FFT window length
    //       Fs:   sample rate in Hz
    //       M:    offset between FFT windows in samples
    //       w:    windowing function; accepted values are
    //             "boxcar", "hann", "hamming", "blackman", 
    //             "bharris", "flattop" and a custom vector
    //             containing N components of the window
    //       grph: graphic output; accepted values are
    //             "lin", "lincmap", "log", "logcmap",
    //             lin/log refers to linear or logarithmic
    //             distribution of magnitude along the
    //             colormap and cmap to the inclusion or 
    //             not of a colormap. The default is no
    //             graph.   
    //       z:    matrix with N/2 rows and
    //             floor((length(x) - N)/M + 1) columns, 
    //             where each column contains the absolute
    //             value of the spectrum beginning at each 
    //             multiple of M
    //       t:    vector containing starting times in s of
    //             each FFT window
    //       f:    vector containing the values of the
    //             FFT frequency bins in Hz
    // 
    // -----------------------------------------------
    // Author: Federico Miyara
    // Date:   2021-07-10
    
    // Argument handling
    rhs = argn(2);
    if rhs<6
        grph = [];
    end
    if rhs<5
        w = "hann";
    end
    if rhs<4
        M = 512;
    end
    if rhs<3
        Fs = 44100;
    end
    if rhs<2
        N = 1024;
    end
        
    // Data for testing purpooses --DELETE--
    if 1==2
        Fs = 44100
        N = 1024
        M = 256
        w = "hann"
        grph = "logcmap"
        tt = 0:1/Fs:2;
        fo = 8000
        fm = 0.5
        df = 4000
        ff = fo + df*sin(2*%pi*fm*tt);
        phi = 2*%pi*cumsum(f)/Fs;
        caso = 2
        select caso
        case 1
            x = 0.4*sin(phi) + 0.3*sin(2*phi);
        case 2
            x = sin(phi) + sin(2*phi);  
            x = min(x, ones(x));          
        end
        wavwrite(x,"d:\hola\fmod.wav")
    end
    
    // Ensure the signal is a column vector
    x = x(:);

    // Number of FFT windows 
    Q = floor((length(x) - N)/M + 1);
    
    // Initialize z
    z = zeros(N/2, Q);
    
    // Generate a matrix whose columns are
    // segments with length N of the signal
    // starting at sample indices multiple 
    // of M
    y = zeros(N, Q);
    for k=1:Q
        y(:,k) = x((k-1)*M+1:(k-1)*M+N);
    end
    
    // Compute window function
    if type(w)==1
        w = w(:);
    else
        select w
        case "boxcar"
            w = 1;
        case "hann"
            w = 0.5 - 0.5*cos(2*%pi*[0:N-1]/N)';
        case "hamming"
            w = 0.54 - 0.46*cos(2*%pi*[0:N-1]/N)';
        case "blackman"
            w = 0.42 - 0.5*cos(2*%pi*[0:N-1]/N)' + ...
            0.08*cos(4*%pi*[0:N-1]/N)'; 
        case "bharris"
            w = 0.35875 - 0.48829*cos(2*%pi*[0:N-1]/N)' + ... 
            0.14128*cos(4*%pi*[0:N-1]/N)' - ...
            0.01168*cos(6*%pi*[0:N-1]/N)';
        case "flattop"
            wft = 1 - 1.90796*cos(2*%pi*[0:N-1]/N)' + ...
            1.07349*cos(4*%pi*[0:N-1]/N)' - ...
            0.18199*cos(6*%pi*[0:N-1]/N)';
        end
    end
    // Replicate window
    ww = repmat(w, 1, Q);
    
    // FFT of the windowed signal along columns
    // and conversion to unilateral spectrum
    z = abs(fft(y.*ww, -1, 1))(1:N/2,:)*2/N;
    z(1,:) = z(1,:)/2;

    // Starting times of the FFT windows 
    t = (0:Q-1)*M/Fs;
    
    // Frequencies of the FFT bins
    f = (0:N/2-1)*Fs/N;
    
    // Graphic output
    cmap = jetcolormap(64);    
    select grph
    case "lin"
        figure();
        clf();
        gcf().color_map= [cmap; [1 1 1]];
        gcf().background = 65;
        Sgrayplot(t, f, z');
    case "lincmap"
        figure();
        clf();
        gcf().color_map= [cmap; [1 1 1]];
        gcf().background = 65;
        Sgrayplot(t, f, z');
        colorbar   
    case "log"
        figure()
        clf();
        gcf().color_map= [cmap; [1 1 1]];
        gcf().background = 65;
        Sgrayplot(t, f, 20*log10(z'));
    case "logcmap"
        figure()
        clf();
        gcf().color_map= [cmap; [1 1 1]];
        gcf().background = 65;
        Sgrayplot(t, f, 20*log10(z'));
        colorbar   
    end
    
endfunction

-------------- next part --------------

// Testing the function spectrogram()

// 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 = 2
select caso
case 1
    // Two harmonics with aliasing 
    x = 0.4*sin(phi) + 0.3*sin(2*phi);
case 2
    // Two harmonics with aliasing and distortion
    // by truncation
    x = sin(phi) + sin(2*phi);  
    x = min(x, ones(x));          
end

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

exec spectrogram.sci

z = spectrogram(x, N, Fs, M, w, grph);