[Scilab-users] spectrogram
Federico Miyara
fmiyara at fceia.unr.edu.ar
Tue Jul 13 06:02:11 CEST 2021
The time-frequency representation implicit in the spectrogram has many
uses apart from sound analysis, for instance ECG, EMG and EEG signals,
mechanical vibration, light signals, seismic signals, sun spots,
circadian rhythms, meteorology, etc.
Here I'm attaching a new version of the function spectrogram() that
includes the possibility of choosing the color map, as in the updated
version of mapsound() proposed by Samuel.
This function allows to choose directly the size of the FFT, the
overlapping (in samples rather than in seconds), the sample rate, the
type of windowing function, the type of graphic presentation, with
linear or log distribution of the levels over the color map, either
showing or not the reference color map, and the color map, either one of
the standard ones introduced as a string or a custom one introduced as a
three-column matrix. The output includes the numerical version of the
spectrogram and, optionally, the time axis and frequency axis, allowing
plotting it in a 3D version with plot3d() as well as postprocessing it
(such as pitch or formant extraction).
I have made no attempt to include other spectrum algorithms apart from
the FFT, or reassignment strategies, which could be a future improvement.
Regards,
Federico Miyara
On 12/07/2021 06:58, Claus Futtrup wrote:
> Hi all
>
> An Amplitude(time,frequency) plot should be named Spectrogram. That's
> the official name for it.
>
> >mapsound(), which plots a spectrogram without numerical information
>
> I didn't check, but if this is true, then an improvement as proposed
> by Federico would be good.
>
> Best regards,
> Claus
>
> On 12-07-2021 11:38, Samuel Gougeon wrote:
>> Hello Federico,
>>
>> As reported at http://bugzilla.scilab.org/16530 , mapsound() could
>> not really be used.
>> It has been rewritten in Scilab 6.1.1.
>> It's help page in PDF is available at
>> http://bugzilla.scilab.org/attachment.cgi?id=5178
>> To be compared to the last 6.1.0 former version
>> https://help.scilab.org/docs/6.1.0/fr_FR/mapsound.html
>>
>> It could be still improved later with additional suggested features.
>>
>> Regards
>>
>> Samuel
>>
>> PS : A redirection mechanism to be set in user's preferences has been
>> implemented in Scilab 6.1.1
>> to automatically redirect help queries in the help browser toward the
>> most relevant page, as uman() does.
>> It has been designed first to help scilab newbies coming from other
>> languages, but can also be used
>> for other common terms.
>> Even if it is not "standard", mapsound is explicit enough, in the
>> right help section, and OK to me.
>> So we could redirect spectrogram to mapsound.
>>
>> Le 12/07/2021 à 01:21, Federico Miyara a écrit :
>>>
>>> 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
>>
>>
>> _______________________________________________
>> users mailing list
>> users at lists.scilab.org
>> http://lists.scilab.org/mailman/listinfo/users
>
>
>
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users
--
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/20210713/5e99e5c7/attachment.htm>
-------------- next part --------------
function [z, t, f] = spectrogram(x, N, M, Fs, w, grph, cmap)
// Spectrogram of a signal
//
// Usage
// [z, t, f] = spectrogram(x, N, Fs, M, w, grph)
// where
// x: signal vector
// N: FFT window length. Default: 1024
// M: offset between FFT windows in samples
// Default: 512
// Fs: sample rate in Hz. Default: 44100
// w: windowing function; accepted values are
// "boxcar", "hann", "hamming", "blackman",
// "bharris", "flattop" and a custom vector
// containing N components of the window.
// Default: "hann"
// 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.
// cmap: a string indicating any of the accepted
// color maps (such as "jetcolormap") or
// any matrix with 3 columns containing the
// RGB description of the colors. Default:
// "jetcolormap"
// 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<7
cmap = jetcolormap(64);
else
if type(cmap)==10
// Accepted color maps
cmaps = ["autumncolormap" "bonecolormap" ...
"coolcolormap" "coppercolormap" ...
"graycolormap" "hotcolormap" ...
"hsvcolormap" "jetcolormap" ...
"oceancolormap" "parulacolormap" ...
"pinkcolormap" "rainbowcolormap" ...
"springcolormap" "summercolormap" ...
"whitecolormap" "wintercolormap"]
if or(cmap==cmaps)
// Get cmap from accepted color maps
execstr("cmap = " + cmap + "(64)");
else
cmap = jetcolormap(64);
end
elseif type(cmap)==1
if isreal(cmap)
if size(cmap)(2)==3 & max(cmap)<=1 & min(cmap)>=0
else
cmap = jetcolormap(64);
end
else
cmap = jetcolormap(64);;
end
else
cmap = jetcolormap(64);;
end
end
if rhs<6
grph = "log";
end
if rhs<5
w = "hann";
end
if rhs<4
Fs = 44100;
end
if rhs<3
M = 512;
end
if rhs<2
N = 1024;
end
// Data for testing purpooses --DELETE--
if 1==2
N = 1024
M = 256
Fs = 44100
w = "hann"
grph = "logcmap"
cmap = "jetcolormap"
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
// Two harmonics with aliasing
x = 0.4*sin(phi) + 0.3*sin(2*phi);
case 2
// Two harmonics with aliasing
// and distortion
x = sin(phi) + sin(2*phi);
x = min(x, ones(x));
end
wavwrite(x,"d:\hello\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;
// Window attenuation compensation
W = sqrt(mean(w.^2));
z(2:$,:) = z(2:$,:)/W;
z(1,:) = z(1,:)/mean(w);
// 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
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 = 3
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));
case 3
Fs = 22050;
t = 0:1/Fs:1*(1-%eps);
x = 0.4 + [1.5*sin(2*%pi*800*t) sin(2*%pi*1200*t)];
grph = "lincmap";
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, M, Fs, w, grph, "jetcolormap");
More information about the users
mailing list