Spectrogram

Obviously neurophysiological signals do not expand in time infinitely. As a matter of fact, the study of their changes (e.g. the detection of the states of consciousness or the electromyographic changes in relation to movement patterns) constitutes the information that we want to extract from our signals.

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.sci

From 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))

Figure 12.8: Figure showing a contour plot of the spectrogram of an electrocardiogram. The bottom frame is the original signal. The left frame plots the spectrum of the signal. The middle frame (the spectrogram) can be used to determine where the different frequencies are located

We can try to make clear the process by calculating the Fourier Transform at three representative points. At each point we apply a hanning window and then we calculate the Fourier Transform.

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:$));

Figure 12.9: This figure tries to show the effect of windowing on a signal. We can see that windowing allows to focus the analysis on some specific segments of the trace and, consequently, the frequencies contained in this particular space can be stressed. The top frame shows the original signal. The second row of frames emphasizes the P wave, the third one emphasizes the QRS complex and the fourth one emphasizes the T component. Notice that the spectrum varies consequently. The vectors containing the spectrum are some of the columns of the matrix represented in the contour plot of the previous figure.

The result shows that the window focuses the Fourier Transform to some components. We can see that the energy of P and T waves is mostly concentrated under 5 Hz while the energy of QRS complex spans from 5 to 30 Hz. To get this result, the choice of the appropriate length of the window is critical since a wider window could mix the frequencies of the different waves while a shorter window could lack frequency resolution to show the changes.