So, we are interested in the description of the modification of the different frequencies with time. An obvious approach is to evaluate the power spectrum of chunks of signal and to represent the power spectrum against time. This kind of representation is called an spectrogram and can be stored in a matrix. Its characteristics are determined by the window that we apply to the signal.
I do not know a direct function that calculates the spectrogram in Scilab but this information is contained in the function `pspect'. Since Scilab is an open source package, we can modify the code to fill our requirements. We want to create a new function (that will be called `spgrm') that returns a matrix containing the power spectrum of each time.
The code of `pspect' can be found in the file `pspect.sci' in the directory `macros/signal/' inside the directory which contains Scilab. I reproduce here some parts of it
function [sm,cwp]=pspect(sec_step,sec_leng,wtype,x,y,wpar) //<sm,cwp>=pspect(sec_step,sec_leng,wtype,x,y,wpar) //Cross-spectral estimate between x and y if both are given //and auto-spectral estimate of x otherwise. //Spectral estimate obtained using the modified periodogram method // x :data if vector, amount of input data if scalar // y :data if vector, amount of input data if scalar // sec_step :offset of each data window // sec_leng :length of each data window // wtype :window type (re,tr,hm,hn,kr,ch) // wpar :optional window parameters // : for wtype='kr', wpar>0 // : for wtype='ch', 0<wpar(1)<.5, wpar(2)>0 // sm :power spectral estimate in the interval [0,1] // cwp :unspecified Chebyshev window parameter //! //author: C. Bunks date: 14 Sept 1988 // Copyright INRIA ... //get number of sections to be used ... //construct window ... //average periodograms sm=0*w; if maxi(size(x))==1 then, ... else, ind=(1:sec_leng); for k=1:nsecs, indd=ind+(k-1)*sec_step*ones(1:sec_leng); xd=x(indd); xe=w.*(xd-(sum(xd)/sec_leng)*ones(1:sec_leng)); fx=fft(xe,-1); if cross==1 then, yd=y(indd); ye=w.*(yd-(sum(yd)/sec_leng)*ones(1:sec_leng)); fy=fft(ye,-1); sm=sm+fx.*conj(fy); else, sm=sm+real(fx.*conj(fx)); end, end, end, sm=sm/(nsecs*wpower); ...The function applies a window to one or two signals and calculates its power spectrum. Finally all the chunks considered are averaged. Its options are the offset of each window, the length of the window and the type of windowing applied. It acts on two signals. The function output is an array containing the average of the power spectrum of all the chunks.
The code begins by calculating the number of sections and creates a window which will be applied progressively to the trace (`w'). Then it generates a vector of zeros of the same dimensions of `w' where the average spectrum will be stored (`sm'). The body of the function calculates the succesive power spectra and averages them with the command
sm=sm+real(fx.*conj(fx));We are going to modify only two lines of the definition of the function and the output of the function. We want that the new function returns all the intermediate power spectra instead of their average. So we are going to edit only two lines
function [sm,cwp]=pspect(sec_step,sec_leng,wtype,x,y,wpar) // Line 1 ... sm=sm+real(fx.*conj(fx)); // Line 94 ...Once we locate the file, we have to copy it to the directory where we are going to work. We copy the file `pspect.sci' (which is included in `macros/signal') in the working signal. Now, in our directory we copy `pspect.sci' to `spectrgram.sci' and we edit it to introduce the changes
cp pspect.sci spectrgram.sci emacs spectrgram.sciFrom our editor, we change the lines `1' and `94' in the next way
function [sm,cwp]=spectrgram(sec_step,sec_leng,wtype,x,y,wpar) // New line 1 ... sm=[sm ; real(fx.*conj(fx))]; // New line 94 ...and we save the changes. We have a new function `spectrgram' contained in the file `spectrgram.sci'. We have made two changes: to modify the name of the function and to store the power spectrum of each segment.
Now it is time to check the functioning or our new creature. We are going to study the evolution of the frequencies present in an electrocardiographic recording sampled at 100 Hz. First of all we are going to load the function in Scilab
-->stacksize(20000000) -->load signals -->ecg = signals(12)(1:300); -->getf('spectrgram.sci');Now we have added a new function to Scilab (`spectrgram') that can be used with the same inputs of `pspect'. Of course the modification tries to minimize the changes and it only produces correct results when we introduce `x' but not `y'.
Since we are interested in changes of frequency inside each beat, we have to use a short window. We can define the degree of overlapping by using secstep. We decided to use a step of one sample. We are going to use a window of length `40' that corresponds to 0.4 s.
-->N = 40; -->fre = (0:N-1)*(1/(N*0.01)); -->ecgspgrm = spectrgram(1,N,'hn',ecg); -->Pxx = pspect(1,N,'hn',ecg);With this code we calculated the frequency array (to plot x-axis in the graphs) from the length of 40 samples at 0.01 s of sampling frequency. Then we use `spectrgram' to calculate the spectrogram and `pspect' to calculate the power spectrum of `ecg' with the same potions (step 0.01 ms, length of window 0.4 ms, hanning window). To represent the values of a matrix we can use different graphs. In this example we use `contour'
xsetech([0,0,1/4,3/4]) plot2d(Pxx ,fre ) xsetech([1/4,3/4,3/4,1/4],[1,-400,275,200]) h= ecg(N/2:$); plot2d(1:length(h),h,1,'002') xsetech([1/4,0,3/4,3/4]) xset('fpf',' '); levels = [1,5,20,60,100]; contour2d(1:size(ecgspgrm,1),fre,ecgspgrm,levels,ones(1,5))
wind = window('hn',N); pos1 = 53; pos2 = 63; pos3 = 93; segment1 = ecg(pos1:pos1+N-1) .* wind; segment2 = ecg(pos2:pos2+N-1) .* wind; segment3 = ecg(pos3:pos3+N-1) .* wind; xset('window',2); xsetech([0,0,1/2,1/4]) plot(ecg(40:120)) xsetech([0,1/4,1/2,1/4]) plot(segment1) xsetech([1/2,1/4,1/2,1/4]) plot2d(fre,ecgspgrm(pos1,1:$)); xsetech([0,2/4,1/2,1/4]) plot(segment2) xsetech([1/2,2/4,1/2,1/4]) plot2d(fre,ecgspgrm(pos2,1:$)); xsetech([0,3/4,1/2,1/4]) plot(segment3) xsetech([1/2,3/4,1/2,1/4]) plot2d(fre,ecgspgrm(pos3,1:$));