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
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
// 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);
_______________________________________________ users mailing list users@lists.scilab.org http://lists.scilab.org/mailman/listinfo/users