How can I use FFT to get frequency of this chemical oscillator (BZ FKN reaction)?

So, I filmed the reaction with my phone and wrote a script to pull the Red, Green and Blue color data from a couple pixels in each beaker as they stirred along with a time stamp in milliseconds.

Using this code I wrote for AutoHotkey:

#NoEnv  ; Recommended for performance and compatibility with future AutoHotkey releases.
; #Warn  ; Enable warnings to assist with detecting common errors.
SendMode Input  ; Recommended for new scripts due to its superior speed and reliability.
SetWorkingDir %A_ScriptDir%  ; Ensures a consistent starting directory.

#SingleInstance Force

A = 297
B = 462
C = 330
D = 454

StartTime := A_TickCount

SetTitleMatchMode, 2
  ControlSend, , {space}, VLC media player
SetWinDelay, 2
SetBatchLines, -1
Loop
{
ElapsedTime := A_TickCount - StartTime
IfWinExist, 20160114_143257.mp4 - VLC media player
{  
WinActivate
  pixelgetcolor, vthiscolor, %A%, %B%, rgb
  vblue := (vthiscolor & 0xFF)
  vgreen := ((vthiscolor & 0xFF00) >> 8)
  vred := ((vthiscolor & 0xFF0000) >> 16)
  pixelgetcolor, vthiscolor2, %C%, %D%, rgb
  vblue2 := (vthiscolor2 & 0xFF)
  vgreen2 := ((vthiscolor2 & 0xFF00) >> 8)
  vred2 := ((vthiscolor2 & 0xFF0000) >> 16)
}

IfWinExist, Pixel 1.txt - Notepad
{  
  WinActivate
}
  sendinput %vred%, %vgreen%, %vblue%, %ElapsedTime%, {Space} {Enter}
  

IfWinExist, Pixel 2.txt - Notepad
{  
  WinActivate
}
  sendinput %vred2%, %vgreen2%, %vblue2%, %ElapsedTime%, {Space} {Enter}
  
}
Esc::
SetTitleMatchMode, 2
ControlSend, , {space}, VLC media player
ExitApp

My script pulls and measures the time the pixels were pulled based on the tick counter of the PC which is time in ms since the PC is started. So it measures time elapsed with

(Start time) - (Ticks since started)

Each pixels rgb value can be r 0-255 g 0-255 and or b 0-255. I've decided to ignore green since the color oscillates between red and blue. So I've plotted r/(r+b) over time to get a %red value.

A picture of the data is attached.

The highest values are when the solution is the most red, and lowest are when the solution is blue. I want information about how often it cycles. The time between changing from red to blue. It's my understanding that this can be done with FFT. I can easily measure the time between cycles just by analyzing the data and getting the time stamps from when my red value is at minimums. However, I'm interested in how the FFT function works and how/if I can apply it to my work.

I am trying to follow along with this guide here: http://www.mathworks.com/matlabcentral/answers/45844-how-to-use-fft-in-matlab-using-imported-data-in-time-domain-excel-file

Here is what I've done/THINK I've done so far:

Import my data as a Nx2 Matrix

Rename my variable to "X" so I can continue using the above guide as a walkthrough.

Note that my Value is 30732 for the variable X, so I will be using the part of the code for an even length:

     xdft = fft(X(:,2));
     % sampling interval -- assuming equal sampling
     DT = X(2,1)-X(1,1);
     % sampling frequency
     Fs = 1/DT;
     DF = Fs/size(X,1);
     freq = 0:DF:Fs/2;
     xdft = xdft(1:length(xdft)/2+1);
     plot(freq,abs(xdft))

Enter Code:

And :(

I've attached my excel sheet.

Thanks, Ryan

Answer 1

[Star Strider]

I cannot understand the file you posted. What is the time and what is the signal?

I used this code, so correct the time ‘t’ and signal ‘r’ references to your file if I have them wrong:

[d,s,r] = xlsread('Ryzx7 Run 4 FFT.xlsx');

t = d(:,3);                 % Time
r = d(:,4);                 % Signal
L = length(r);

Ts = mean(diff(t));
Fs = 1/Ts;
Fn = Fs/2;

Fr = fft(r)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:length(Fv);

figure(1)
loglog(Fv, abs(Fr(Iv)))
grid
axis([0  180    ylim])
xlabel('Frequency (Hz)')
ylabel('Amplitude')

[Ryzx7]

The excel workbook contains 3 worksheets (the tabs on the bottom)

Sheet labelled Pixel 1 HD:

Each row is a measurement of the selected pixel at time (T)

Columns A B C labelled (R), (G) and (B) are the color values in decimal form for Red Green and Blue of that pixel at that given time in milliseconds (T). The signal is labelled (%R) and is in column E of the sheet. The (%R) column is the red value of a given pixel divided by the sum of the red and blue values of that given pixel so it's value is part red divided by parts red+blue or R/(R+B). It's the measure of how red the solution is.

Pixel 2 HD: is the same as Pixel 1 HD except I'm sampling a different pixel from the same beaker.

The sheet labelled T %R is the time column from Pixel 1, and the %R columns from pixel 1 and pixel 2. The %R columns are averaged to try and get a better idea of the color of the whole solution. So on the 3rd sheet labelled "T %R", I have Time (T) the average redness of my solution from 2 pixels. The final column was just the difference in time because my script doesn't sample pixels at exactly 1 sample per.125 seconds.

Thank you for your code, I used it to generate this:

I'm not sure what it means, but it's much more interesting than what I had. I thought it would like like a flat line with a spike at ~33 Hz corresponding to the way the solution turns blue every 33 seconds. Something like this kind of action:

it is my understanding that

[d,s,r] = xlsread('Ryzx7 Run 4 FFT.xlsx')

means create 3 variables d: A matrix without headers s: the text based headers r: A cell that contains text and numbers

so it's just a fast way to digest the excel sheet into its parts but we really only need d for the FFT. This command pulls from sheet 1 of the workbook by default I guess.

t = d(:,3);

Says: create a variable t that is equal to the contents of the matlab variable d which is the matrix we imported from excel and inside that is an ordered pair. The colon means ... I don't know what the colon means. The 4 means the fourth column which in this case is my time variable. So I change the 3 to a 4. The variable r has the same format, but r the signal is going to be %R and that's the 5th column so I change that to a 5.

Ts: that's like my dT column from sheet 3. So that's taking the average of the differences of each T to get the average change in T?

Fs: Frequency of sampling so 1/Time between samples

Fn: Fn something to do with the FFT needing a frequency that's half the frequency sampling

Fr the FFT applied to r (my signal) divided by the number of values in the r column (L, length of r)

Fv: something to do with this that I don't understand yet.

Iv: same as above

Create a figure

log(log(Fv, abs(Fr(Iv))) grid: axis limits not sure on this syntax or the 180 limit labels

[Star Strider]

I couldn’t open your file in Excel (that is a problem with Firefox), so I had no idea what it contains. That’s the reason I asked. I opened this thread in IE and successfully opened the spreadsheet in Excel.

Now that I know it has three sheets and the third is the one I need to look at, I changed my code accordingly. I always ask for all three outputs from xlsread so I know what the structure of the relevant sheet is, and what column is what, so I don’t have to open it in Excel as well.

Your data are not absolutely regularly sampled, but are close enough. Ideally, I might use the Signal Processing Toolbox resample function, or possibly interp1 to get a regular sampling frequency, but yours are close enough to regularly-sampled. The mean difference between samples is 0.1246 and the standard deviation is 0.0137. (The code I used for the fft analysis is from the R2015a documentation that is easier to understand: fft.)

‘Fr’ is the fft of your signal (‘%R’). It is normalised for length and doubled because we are looking at only half of the fft, so the amplitudes are correct if we multiply it by 2.

‘Fn’ is the Nyquist frequency. It is the highest resolvable frequency in a sampled signal.

‘Fv’ is the frequency vector corresponding to the fft up to the Nyquist Frequency.

‘Iv’ is the index vector, so that the plotted ‘Fr’ will match the length of ‘Fv’.

The new code is:

[d,s,r] = xlsread('Ryzx7 Run 4 FFT.xlsx',3);    % Read Third Sheet Of File

t = d(:,1);                                     % Time Vector
r = d(:,2);                                     % Data Vector
L = length(r);

Ts = mean(diff(t));                             % Sampling Interval
Fs = 1/Ts;                                      % Sampling Frequency
Fn = Fs/2;                                      % Nyquist Frequency

Fr = fft(r)*2/L;                                  % FFT Of Signal (Normalised For Length And Amplitude)
Fv = linspace(0, 1, fix(L/2)+1)*Fn;             % Frequency Vector
Iv = 1:length(Fv);                              % Index Vector

figure(1)
semilogy(Fv, abs(Fr(Iv)))
grid
axis([0  Fn    ylim])
xlabel('Frequency (Hz)')
ylabel('Amplitude')

The plot gives peaks at about 0.8, 1.0, 1.6, and 2.0 Hz. It’s appropriate to consider the last two as harmonics of the first two. You can just use plot instead of semilogy, but here it was easier to see the peaks with semilogy.

Note that loglog is not log(log(x)). The loglog plot uses log axes on both x and y making low frequencies easier to resolve. That is not necessary here, so I changed the code accordingly.

[Ryzx7]

Star Strider, thank you so much for taking the time. I ended up with this as my output:

Had a little bit of a struggle with the xlsread, but with your code I was able to easily import my data into a Nx2 matrix and proceed from there. Really appreciate you taking the time to help me with this. Thank you.

[Star Strider]

As always, my pleasure.