Image analysis of brain image

I have a small project to do on Image Analysis. I have an excel file for different 6 different pixels and frames for intensities of a brain signal. My job is to analyze the theta waves which in 2x3 pixels region. These are my following task:

1. to compute fft
2. to find frequency for theta waves
3. to use filter to remove noise

So what I think is for 1, I need to create a 2d matrix to load the data for pixel 2 and 3, and compute fft2(). However, I have only done a couple of examples for 1D fft yet. I am not even sure if we need to use 2d fft for theta waves analysis. I am completely stuck. Can somebody help me getting started?

Answer 1

[William Rose]

@reema shrestha

I assume that you have a record of intensity (or voltage, if it is an EEG) as a function of time for each pixel. In other words, you have the equivalent of a gray-scale video with 2x3 pixels. It woudl help if you upload your Excel file so readers understand the problem better. What s the smapling rate and how manysamles or seconds of data do you have?

If I am right that you have a record of amplitude versus time for each pixel, then you can do six 1D FFTs: one for each pixel. You say you have to find the frequency of theta waves. Theta waves are in the 4-8 Hz band. Once you hav computed the fft, you have several options for esitmating the theta wave frequency, including but not limited to the following

• Find the frequency in the 4-8 Hz range that has max power (or, equivalently, max amplitude).
• Take a power-weighted average frequency in the 4-8 Hz band.
• Find the median frewuncy in the 4-8 Hz band.

The above approach, if done separately for each pixel, will probably give a different result for each pixel. If you want to combine information from all the pixels, you could average the theta wave frequencies from each of the pixels to get an average frequency.

An alternative would be to coombine the image data from the six pixels into a single channel or signal. Then you just do one FFT on the one signal, and find the theta wave frequency using one or mor of the methods which I listed above.

As for filtering: you can filter each pixel indidually or you can combine them and filter the combined signal. Combining the pixels is in itself a form of filtering, since it averages 6 into 1. If the noises on pixels are somewhat uncorrelated, then averaging them reduces the noise.

[reema shrestha]

This is a small part of my data. the samples are taken each 77ms. Regarding the 2x3 pixel, there is a small region in the image, for which I need to calculate the frequency of the theta wave in that region.

To be more precise, the instruction is "This assignment will detect the theta waves in the right hippocampus of the brain. A region of 2 x 3 pixels, as pointed by the arrow in the figure above, is used for this purpose.

work:

To produce frequency spectral of the theta waves using FFT. To identify the frequency of the theta waves in the above chosen region." so that means i need to plot two ffts. For the first one, i need to combine all the pixel, and for the second one, i need to combine only 2 and 3 pixel. or do i have to only deal with 2x3 pixel?

[William Rose]

@reema shrestha,

Esitmating the EEG spectrum by doing the FFT is simple, but the amplitude estimate at each frequency is likely to be somewhat noisy. I assume you are aware that, if X is a real vector, fft(X) returns a complex vector which is the two-sided discrete Fourier transform. You have to convert it to an amplitude spectrum or power spectrum and you have to compute the correspoding vector of frequencies. The noisiness of the spectral estimate at each frequency will be lower if you use some kind of signal averaging technique. The Welch spectrogram is the most widely used signal averaging technique for spectra. Matlab has a nice built-in routine for this: pwelch.

https://www.mathworks.com/help/signal/ref/pwelch.html

The p in the name is because pwelch computes an estimate of the power spectrum, rather than the amplitude spectrum. That is fine for your work, since your question implies that you are not interested in phase information. You can use the power spectrum from pwelch to find the theta wave frequency, using one of the methods I suggested in my other posting.

[reema shrestha]

Thank you for your detailed answer.

[William Rose]

@reema shrestha,

The instructions you recieved are not very clear, at least not to me. It sounds like you have a Row by Col by N dataset, and you are to analyzea 3x2xN subset. You did not say what N is (how many frames). My interpretation of the instructions you have quoted is that you are to combine the 3x2 pixels into one. For this, you might as well average them together. This results in a Nx1 vector, which you may call X. If N is not even, I recommend deleting the first or last row, so that N is even. This will simplify the extraction of the 1-sided spectrum from the 2-sided spectrum.

dt=0.077;
Y=fft(X);
A2 = abs(Y/N);
A1 = A2(1:N/2+1);        %we assume N is even
A1(2:end-1) = 2*A1(2:end-1);
f = (1/dt)*(0:(N/2))/N;
plot(f,A1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|A1(f)|')

The code above is mostly taken from an example here. It computes and plots the single-sied amplitude spectrum. Then you find the theta wave frequency in the 4-8 Hz range, using one of the approaches I mentioned earlier.

[reema shrestha]

Following your suggestion, I filtered each 6 pixels first with both FIR amd IIR filter and computed the fft for each pixels. Then later, averaged the amplitudes of the ffts. This gave me the following output. For the frequency of theta wave, I chose to use the max amplitude value, even though the graph is not smooth. But I am still not sure, if it is correct. The peak for theta wave seems small compared to some of the peak near 1 Hz. Also, I have a question. While computing the single spectrum fft, the max range of frequency/Nyquist freq is less than 8 Hz which is the range of theta wave. If we were taking the two sided spectrum, it wouldn't have been the case, I guess. I'd appreciate your opinion for this following output.

Answer 2

[William Rose]

@reema shrestha,

You are right about the limited frequency range. You said the sampling interval was 0.077 second, i.e. sampling rate fs=13.0 Hz. The Nyquist frequency is fNyq=fs/2=6.5 Hz. The amplitudes in the FFT at frequencies from 6.5 Hz to 13.0 Hz (the top half of the two-sided spectrum) are just the "folded over" amplitudes from the 6.5 to 0 Hz range - IF the analog signal does not have any power above 6.5 Hz. The reason this happens is due to the mathematics of discrete sampling and sinusoidal signals. If there is power above the Nyquist frequency, then the FFT at all frequencies in th two-sided spectrum is a complicated overlay of the amplitudes below and above the Nyquist frequency. You can never disentangle them, after they have been sampled. This is called aliasing, because a signal at a high frequency takes on the appearance, or identity, of a signal at a lower frequency.

Therefore you have two problems:

1. The 13 Hz sampling rate is not fast enough to capure the 4-8 Hz band, since the associated Nyquist frequency is 6.5 Hz, which is below 8.
2. Aliasing almost certainly occured when the signal was sampled. I say this because the original amplitude spectrum plot shows a fairly constant amplitude (+- some noise) from 2 Hz up to the Nyquist frequency. It very unlikely that the analog signal was sharply band limited right at 6.5 Hz. It is much more likely that the signal had power above 6.5 Hz and that the amplitude we see below 6.5 Hz is a mix of signal that really was below 6.5 Hz plus analog components above the Nyquist frequency, whch got aliased down into the 0-6.5 Hz range due to sampling at 13 Hz.

If this data were from my advisor, I would tell them that it is not possible to reliably determine the theta wave frequency or theta wave content from this recording, due to the low sampling frequency and the likelihood of aliasing. I would propose that we collect more data at a much higher video sampling rate. (See below.) But we can still try to analyze the spectrum that you have, from 4 to 6.5 Hz.

I agree that the frequency which has maximum amplitude may not be a great estimate of mean theta wave frequency, because the spectrum is quite spiky or noisy. Thefore the frequency of one skinny tall spike may not be a very reliable indicator. I would compute the mean and median frequencies as well.

Sinc eyou have 115.5 seconds of data (1500 x 0.077), I recommend you try pwelch() . You will lose frequency resolution (i.e. in your spectrum will be bigger) but the error bars on your spectrum will get smaller. If you want to post some data, such as the raw data array (3x2x1500), or the orignal one-sided average spectrum (should be a vector of amplitudes, length 751), that could be helpful.

Sampling rate

To avoid aliasing, the analog signal needs to be analog-filtered to remove power above the Nyquist frequency, before it is sampled. The original amplitude spectrum plot which you show makes me suspect that this signal was not sufficiently low-pass filtered before it was sampled.

The sampling rate is chosen based on 1. the range of frequencies you wish to study (up to 8 Hz in this case) and 2. on the properties of the analog low pass filter used to prevent aliasing. This is its own complicated subject, but the end result is that you need to sample fast enough that the Nyquist rate is well above the highest frequency of interest, often by a factor of 2 or more. To study frequencies up to 8 Hz, I want a Nyquist frequency of 16 (or more, depending on my anti-aliasing filter), and therefore a sampling rate of 32 Hz or more. Most EEG system sample at 240 Hz, it says here.

Your signals comes from images. There is no pracitcal way to analog-low-pass-filter the light intensity from a scene. (That is why aliasing causes helicopter blades in a video image to appear to be going very slowly. Their high frequency has been aliased down to a much lower frequency by video sampling.) Therefore it is good to sample at the highest rate which your video system allows.

[reema shrestha]

So, this is my data. It's just a part of an assignment. Thus, there is nothing I can do with the aliased signal. So, according to you, the spikes on the graph is caused due to aliasing, right? About trying pwelsch, one of the objective of this assignment is to compute frequency in this region, and if pwelsch loses the frequency resoultion, is it even worth it to try it? Since, I am stuck with this data, I guess, the only option left for me is to take mean and median of the frequency?

The only thing I am confused right now is, there is a small table given at the right side in the excel for the ordering of pixel. Since we are taking the average of all the six pixels, I guess, there is no difference/use of that table. Right?

[William Rose]

@reema shrestha

The table in the Excel sheet for pixl order is not important. It tells you which column of the spreadsheet corresponds to which pixel.

In Excel, I computed the mean of 6 pixels intensity at each frame, and save the result as a single clumn of 1500 numbers, in a text file. The file is attached.

The attached Matlab script reads the data from the text file, does the FFT, plots the single sided amplitude spectrum, and finds the peak frequency, average frequency, and median frequency within the 4-8 Hz band. The band ends at 6.5 Hz, due to the 13 Hz sampling rate. The mean and median frequencies are 5.23 Hz. This is about halfway from 4.0 to 6.5, which is not suprising, since the amplitudes are fairly constant over the band. If the amplitudes were exactly constant, the median and mean would be exactly at the halfway point of the band. The frequency of max amplitude is 5.65 Hz.

[reema shrestha]

Thank you for the detail answer.