Circadian analysis using FFT()
10 views (last 30 days)
Show older comments
I need to get the power spectrum of gene expression data that were sampled every 4 hours for a 48 hour period (13 samples, including zero point). I did an initial analysis based on a sound recording example but I think I'm missing something.
This is the code I'm using:
signal = load('data.txt')
N = length(signal);
fs = 1/(4*60*60); %Samples per second
fnyquist = fs/2; %Nyquist frequency
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N;
N_2 = ceil(N/2);
plot(fax_Hz(1:N_2), X_mags(1:N_2))
I want to know if there are frequencies associated with the period of the signals. This data should have a period around 24 hours.
Mike
0 Comments
Accepted Answer
Star Strider
on 19 Oct 2014
Edited: Star Strider
on 19 Oct 2014
I would do it slightly differently, and use your 4-hour sampling time directly:
GE = dlmread('Circadian analysis using FFT_data.txt');
T = [0:length(GE)-1]*4; % Sampling Times (Hours)
Ts = 4; % Sampling Interval (Hours)
Fs = 1/Ts; % Sampling Frequency (1/Hour)
Fn = Fs/2; % Nyquist Frequency (1/Hour)
FGE = fft(GE)/length(GE);
Fv = linspace(0,1,length(FGE)/2+1)*Fn;
Iv = 1:length(Fv);
figure(1)
% plot(Fv, abs(FGE(Iv))) % Frequency = 1/Hour
plot(Fv*24, abs(FGE(Iv))) % Frequency = 1/Day
grid
xlabel('Frequency (D^{-1})')
% xlabel('Frequency (H^{-1})')
ylabel('Amplitude')
The commented lines in the plot are in frequencies of 1/Hour. The uncommented lines are in frequencies of 1/Day, or (obviously) 1/(24 Hours).
EDIT to add plot:
6 Comments
Star Strider
on 20 Oct 2014
My pleasure!
Signal processing techniques do not care what sampling times or time base you use so long as you are consistent. If you are acquiring your data in terms of hours rather than seconds, it makes sense to do the analysis in units of hours. (One day is 86400 seconds, so that would create some extremely inconvenient units if you defined your frequencies in Hz.)
You can probably do many calculations in MATLAB that you can do in R, although they may require Toolboxes that are not part of core MATLAB. I will help as much as I can with them if you want to expand your use of MATLAB.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!