Fourier analysis/s​pectrum/FF​T/transfor​m Problem

Hello all.

I am studying a system and want to investigate how the frequency of y(3) varies under different conditions. However, my the fft I perform on it tells me the frequency is zero, which must be incorrect. I have tried a stack of things but can't see what the problem is. I'm relatively new to the matlab environment.

The main code:

close all; clear all; 

%--------------------------------------% 
Time control and initial conditions

timeinmin=600; % Time in minutes

tfin=timeinmin*60; % Time in seconds

%--------------------------------------% 
Initial Levels and system Saturation

totIKK=0.08; % Total IKKn at start

totNFkB=0.08; % Total NFkB cytoplasmic at start

T=1;         % Saturation of the system

yint=[totNFkB,0,0,0,0,0,totIKK,0,0,0,0,totIKK,totNFkB,T]; % Initial conditions 

%--------------------------------------% 
ODEs

[t,y]=ode45(@nfkb,[0 tfin],yint); % calls Model

%--------------------------------------% 
Extract Data from the ODE 

fs = 1; % The sampling freqency in Hz

time = 0:1/fs:tfin; % Define time of interesrt 

Datapoints = tfin*fs; % # of Data points being used

yi = interp1(t,y(:,3),time);

%plot(time/60,yi,'o',t/60,y(:,3))  % test to ensure it is working 

%dbstop 

%--------------------------------------% 
Fourier Analysis from 

NFFT = 2^nextpow2(tfin); % Next power of 2 from length of y

Y = fft(yi,NFFT)/Datapoints; 

f = fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.

plot(f,2*abs(Y(1:NFFT/2+1))) 

title('Single-Sided Amplitude Spectrum of y(t)')

xlabel('Frequency (Hz)')

ylabel('Amplitude - |Y(f)|')

% Get value of frequency

[B,IX] = sort(2*abs(Y(1:NFFT/2+1))); 

BFloor=0.1; %BFloor is the minimum amplitude value (ignore small values) 

Amplitudes=B(B>=BFloor) %find all amplitudes above the BFloor 

Frequencies=f(IX(1+end-numel(Amplitudes):end)) %frequency of the peaks

The ODEs:

function dydt=nfkb(t,y)

parameters;

dydt=[kd1*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-ka1*y(2)*y(1)-ki1*y(1)+...

      kdegc*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))+ke1f*y(3)+kdegpin*y(11);

      kd1*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-ka1*y(2)*y(1)-ki2*y(2)+...

      kdegf*y(2)+ktria*y(6)-kc2*y(8)*y(2);

      kd1*y(5)-ka1*y(4)*y(3)+kv*ki1*y(1)-kv*ke1f*y(3);

      kd1*y(5)-ka1*y(4)*y(3)+kv*ki2*y(2)-kv*ke2*y(4)-kdegf*y(4);

      ka1*y(4)*y(3)-kd1*y(5)-kv*ke1c*y(5);

      kitria*(y(3).^h)/((y(3).^h)+(k.^h))-kdegtia*y(6);

      kp*(y(12)-y(8)-y(7))*kbA20/(kbA20+y(10)*y(14))-y(14)*ka*y(7);

      y(14)*ka*y(7)-ki*y(8);

      kitra*(y(3).^h)/((y(3).^h)+(k.^h))-kdegta*y(9);

      ktra*y(9)-kda*y(10);

      kc3*y(8)*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-kdegpin*y(11);

      0;

      0;

      0];

The parameters:

%rates

kv=3.3; % C:N ratio

kp=0.0006; % Conversion rate of IKKi

ka=0.004; %activation rate of IKK

ki=0.003; % Usage rate of IKK

kda=0.0045; % Deg rate of A20 protein

kbA20=0.0018; % A20 concentration

kc2=0.074; % phosopho of IkB

kc3=0.37; % phospho of IkB-NFkB

kdegpin=0.1; % splitting of Ikb-NFkB due to phospho

ka1=0.5; % nuclear IkB-NFkB formation 

kd1=0.0005; % dissisociation of nuclear IkB-NFkB 

kdegf=0.0005; % Deg rate of nuclear IkB

kdegc=0.000022; % dissisociation of cyto IkB-NFkB

ki1=0.0026; % nuclear entry of NF-kB

ke1f=52*10^(-6); % nuclear exit of NF-kB

ke1c=0.01; % Transport of IkB-NFkB out of nuc.

ki2=0.00067; % nuclear entry of IkB

ke2=3.35*10^(-4); % Nuclear exit of IkB

h=2; % hill coefficient of IkB transcription

k=0.065; % dissasociation constant for IkB tanscription with NFkB

kitria=1.4*10^(-7); % trans rate of IkB

ktria=0.5; % translation of IkB

kdegtia=0.0003; % Deg rate of IkB RNA

kitra=1.4*10^(-7); % Transcription rate of A20 RNA

ktra=0.5; % transcription rate of A20

kdegta=0.00048; % deg rate of A20 RNA

Thank you in advance for any insights.

[Dr. Seis]

What do you mean the frequency is zero? When you run this, what is/are the value(s) associated with "Amplitudes" and "Frequencies"?

[Calum]

Let's say we run the code with the initial conditions stated above. An amplitude of 0.1048 and frequency of zero is found. However, I feel the problem lies firstly in the spectral analysis. If one plots y(3) as a function of time it is clearly an oscillatory system, but the the spectral analysis I perform tells me there are no frequencies present. Thus something I can't understand is going wrong in my approach at the FFT step.

Answer 1

[Rick Rosson]

Please format your code. If you don't know how, click on "Markup Help".

[Calum]

Hope that is what you meant.