Frequency analysis of circadian data

3 ビュー (過去 30 日間)
Gabriela Da Silva André
Gabriela Da Silva André 2019 年 2 月 21 日
編集済み: Gabriela Da Silva André 2019 年 2 月 27 日
Hi everyone
i'm analysing the body core temperature, which should follow the circadian rythm. With a frequency analysis i want to find out the period of my data, to be able to fit a cosinor over my body core temperature curve. I have used adapted a code from a user in here, but i don't get the results i would expect (period should ± around 24h). The problem also is, that whenever i adapt my data ( at the moment i have data every minute), for example adapt it to data every 5 minutes, i get another result, which is odd because i would expect a similar result.
%% Read data
if isempty(path)==0
CoreTempPath=uigetdir('/Users/gabriela./Desktop/');
end
CoreTempData=CoreTemp_Open(CoreTempPath);
[l,j]=size(CoreTempData);
%% Applied Filters
% hampel filter: median filter
for i=1:j;
[TempFilt,k]=hampel(CoreTempData(i).temp,50);
CoreTempData(i).k=k;
CoreTempData(i).tempFilt=TempFilt;
end
% moving average filter
a=1;
windowlength=25;
b=1/windowlength.*ones(1,windowlength);
for i=1:j;
TempFilt2=filtfilt(b,a,CoreTempData(i).tempFilt);
CoreTempData(i).tempFiltAverage=TempFilt2;
end
This are the first steps i do to read the data in. (there are for loops because i use to read several body core temperature data at once.
Then i average my data, so i get a data point every minute, with the following code
for i=1:j;
tt=timetable(CoreTempData(i).dateTime,CoreTempData(i).tempFiltAverage);
tt2=retime(tt,'minutely','mean') % mean every minute
CoreTempData(i).meantempwm=tt2.Var1; %mean Temperature over 1 minute with missing data
CoreTempData(i).meantime=tt2.Time;
timetemp=fillmissing(tt2,'linear');%interpolating missing data
CoreTempData(i).meantemp1=timetemp.Var1;
end
then i used the following code to do the frequency analysis and find the period
%% Frequency Analysis using fft --> power vs frequency
for i=1:j;
n = length(CoreTempData(i).meantemp1);
y = fft(CoreTempData(i).meantemp1);
y(1) = []; %sum of the data and can be removed
%code addapted. source: https://ch.mathworks.com/help/matlab/examples/using-fft.html
n = length(y);
power = abs(y(1:floor(n/2))).^2; % power of first half of transform data
nyquist = 1/2; % nyquist
freq = (1:n/2)/(n/2)*nyquist; % equally spaced frequency grid
figure
plot(freq,power)
xlabel('Cycles/minute')
ylabel('Power')
figure
period = (1./(freq));
plot(period,power);
xlabel('min/Cycle')
ylabel('Power')
hold on;
index=find(power==max(power));
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25);
text(period(index)+2,power(index),['Period = ',mainPeriodStr,'min']);
hold off;
end
The result i get isn't quiet that what i would expect. But the thing is i don't know where the mistake in my code is. So if anyone could help, i would be very grateful!
Thank you

回答 (0 件)

カテゴリ

Help Center および File ExchangeFourier Analysis and Filtering についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by