Digital signal processing - average fit of periodic signal?

Dear MATLAB Central,

I would like to process periodic signals using Matlab. I uploaded an example csv file with several periodes: http://idisk.mpi-cbg.de/~weber/matlab/e005_ortho_plane116_roi1.csv

Here are three of the periodes plotted:

My goal is to get an average fit for all the periodes, the overlay of all curves - the "ideal" curve. What is your opinion on the best strategy? I first need to get the frequency (which could change to a small extend), then move the curves above each other using the calculated frequency, and average them.

Any input appreciated! :)

Best regards, Michael

P.S.: I am using MATLAB 2010a with the image processing tool box.

Answer 1

[Wayne King]

I would fit a linear model in the frequency domain. You have one clearly dominant periodicity in your data with a frequency of 17 cycles per 500 samples. I don't know what your sampling interval is, so I can't say anything more specific about the period. To get a better fit overall, I'll use the 3 largest periodicities in your data, but as I said, one oscillation is clearly dominant.

In my example, I have labeled your data as planet_data. So rename your data planet_data in your workspace, just the numeric vector, and you should be able to copy and paste the below.

 ts = detrend(planet_data,0);
 tsdft = fft(ts);
 freq = 0:1/length(ts):1/2;
 plot(freq,abs(tsdft(1:length(ts)/2+1)),'bo-','markerfacecolor',[0 0 1]);
 xlabel('Cycles/Sample'); ylabel('Magnitude');
 title('Discrete Fourier Transform Magnitudes');
 mu = mean(planet_data);
 tsfitdft = zeros(500,1);
 [Y,I] = sort(abs(tsdft),'descend');
 indices = I(1:6);
 tsfitdft(indices) = tsdft(indices);
 tsfit = ifft(tsfitdft,'symmetric');
 tsfit = mu+tsfit;
 figure;
 plot(planet_data); hold on;
 plot(tsfit,'r','linewidth',2)
 title('Time-domain Fit');
 xlabel('Sample'); ylabel('Y-value');

[Michael]

Dear Wayne,

thanks a lot for your input, works very well so far! What I am not completely confident with is the final fit. The inverse FFT of 'tsfitdft' results in a periodic function with varying amplitude. Could that be because of the way you compute the three largest periodicities? They are 17, 18 and 35 but should rather be the three distinct peaks at 18, 35 and 51. I will play around with the code. Would you say that an increased sampling rate and/or sampling length would help to get better results? I read that the sampling length has a high impact on the results of the FFT.

[Wayne King]

The time-varying amplitude that you see is likely just the result of the sine waves interacting (the phase effects). If you just fit the mean plus the largest Fourier coefficient (in magnitude), that will not vary in time.

It is quite possible if you pad the DFT, you will interpolate the peaks to land directly on a DFT bin. Alot depends on how much you know about what periodicities you expect. You never said what your sampling rate is, or what knowledge you have about which periodicities might be present, if any.

[Michael]

The data sets are the result of calcium measurements inside zebrafish hearts. Those measurements are carried out using fluorescence microscopy of a calcium-sensitive marker. Calcium itself induces the contraction of the heart muscle, and those contractions are periodic. The presented curve shows the intensity changes in one particular region of the heart wall. What I would like to obtain from those measurements are the distinct shapes of the curve, and ultimately the pattern of the conduction wave as it travels along the heart wall. The sampling rate for the example data set is 63.2 Hz (the frame rate of the camera). I will soon use another camera with a much higher frame rate (200-400 Hz), and can possibly acquire longer movies. The frequency of the conduction should be rather stable over a certain time, as the conduction is controlled by rhythmic pacemaker cells (and quite crucial for a constant blood flow).

Answer 2

[Jürgen]

if you do not have the curve fitting toolbox then your approach sounds exceptable to find the average signal , put all values for a periode in row in a matrix and calculate the average over the columns, regards,Jürgen

[Michael]

Dear Jürgen,

I read about the curve fitting toolbox before. Are you yourself using it? My understanding is that it mainly helps to manually find the best fit - a bit like trial & error?

My problem with simply averaging all the curves is that I would most probably 'wash out' the distinct shape, and that would be a pity. I guess I have to interpolate each individual curve before doing the averaging.

[Jürgen]

I would not say that it is manually or trial and error, but yes it provides you with tools the look at differents fit. Actually since I have the toolsbox I do not know exactly with function are in the basic matlab and which are in the toolbox regardsJ

Answer 3

[Star Strider]

A File Exchange function that could help you is Peak Fitter. It could help you identify the locations and other characteristics of the peaks. With that information, you could then overlay them and average them. (I have not yet used it myself.)

[Michael]

Thanks for the tip Star Strider, I will give the Peak Fitter a try!

[Star Strider]

After learning that the curve was zebrafish heart calcium concentrations, a brief search of the available literature revealed that the waveform is essentially that of a relaxation oscillator. I could not find an explicit method in the papers I could access, so I devised one of my own. (Wavelets are likely a more robust solution, especially for large data-sets where you may also want to compare waveform characteristics between experiments.) This is a condensed version of my original code, but it should work. The strategy is straightforward, if kludgy:

• Define an acceptably representative exponential function,
• Select a representative section of the data to fit it to (I used lsqcurvefit), and once fit and identified,
• Convolve it with the data (for some reason the built-in convolution functions did not give acceptable results, so I created my own version),
• Identify the ends of each beat,
• Create an ensemble matrix from the original data vector,
• Calculate descriptive statistics

The most difficult part of the code was making the legend in figure(7) work the way I wanted it to. That's not critical, but I want the plot to look decent. (The documentation is confusing.)

The filled contour plot in figure(8) provides some confirmation that the routine works as I want it to.

I created a ‘.png’ of figure(7), but since there is no way to upload it directly to ‘Answers’ without using some third-party upload site (I looked through all the documentation and previous ‘Answers’ posts on that topic I could find), I decided not to post it.

I do not pretend that my code is the best or only solution to your problem. I can only claim that it seems to provide what in my experience appear to be reasonable results.

% ==========================  START CODE  ==========================

Grafont = 'Calibri';

% READ DATA FILE: 
fidi = fopen('e005_ortho_plane116_roi1.csv')
Data = textscan(fidi, '%*3d e005_orthogonal_plane116-1.tif: %f %f', 'Delimiter',',', 'HeaderLines',1 )

Time0 = cell2mat(Data(:,1));
CaConc0 = cell2mat(Data(:,2));
if CaConc0 > 3200
    CaConc0 = 3200;
end

Fs = 63.2;          % Sampling frequency (Hz)
Ts = 1/Fs;          % Sampling time (sec)

% GENERATE A PROTOTYPICAL EXPONENTIAL FUNCTION: 
Ts1 = [0:32]';
fexp = @(p,x) p(1)*sin(2*pi*p(2)*x).*exp(-p(3)*x);
Parameters = [  1.056503647716452E+03;    1.538886399814761E-02;    8.267074622465402E-02];
% (Plug ‘Parameters’ into ‘fexp’ to create ‘fitxp’)
fitxp = fexp(Parameters,Ts1);

% DETREND AND AUGMENT THE ‘CaConc0’ VECTOR TO CREATE CaConc1: 
CaConc1 = detrend(CaConc0);
[MinCaConc0, idxmin] = min(CaConc0);
CaConc1 = CaConc1+abs(min(CaConc1))+1;          % Be sure CaConc1 > 0

CaConc1 = [CaConc1;  zeros(fix(length(fitxp)/4),1)];    % AUGMENT TO SET UP FOR KLUDGY CONVOLUTION

% DO KLUDGY CONVOLUTION (TO IDENTIFY DEPOLARIZATIONS AND DETERMINE PEAKS
%   AND VALLEYS)
nmfitxp = norm(fitxp);                          % 2-Norm of ‘fitxp’
wndw = [1:length(fitxp)]';                      % Define index range for convolution
for k1 = 1 : (length(CaConc1)-length(wndw))
    CaConcIdx = wndw+(k1-1);                    % Increment window by 1 for each iteration
    nmCaConc1 = norm(CaConc1(CaConcIdx));       % 2-Norm of ‘CaConc1’ segment
    dp(k1) = dot(CaConc1(CaConcIdx),fitxp);     % ‘dp’ is dot product of ‘fitxp’ and an equal-length segment of the incremented CaConc1 vector 
    cc(k1) = dp(k1)/(nmCaConc1*nmfitxp);        % ‘cc’ is an ersatz correlation coefficient (for later computational convenience) 
end

% FIND PEAKS AND VALLEY AND DETERMINE THE MAXIMUM DISTANCE BETWEEN THE
%   VALLEYS
[Pk,PkIdx] = findpeaks(cc, 'MinPeakDistance',20);               % IDENTIFY PEAKS
[Vl,VlIdx] = findpeaks(max(cc)+1-cc, 'MinPeakDistance',20);     % IDENTIFY VALLEYS

DifVlIdx = diff([0; VlIdx']);   % Find intervals between VALLEYS (Start Points for each contraction)

% CREATE AN ENSEMBLE MATRIX OF THE IDENTIFIED CALCIUM FLUX EPISODES: 
CaConcEnsb = zeros(max(DifVlIdx),size(VlIdx,2)-1);      % Preallocate
for k1 = 1:length(DifVlIdx)-1
    IdxRng = VlIdx(k1):VlIdx(k1+1);                     % Define Index Range From ‘Valley’ Indices
    LenColIdxRng = length(IdxRng);                      % Calculate Length of Each Index Range (Varies by Depolarization) 
    CaConcEnsb(1:LenColIdxRng, k1) = CaConc(IdxRng);    % Create Next Ensemble Entry
end

% CALCULATE STATISTICS ON THE ENSEMBLE: 
CaConcEnsbMean = mean(CaConcEnsb,2);
CaCondErr95 = [CaConcEnsbMean-1.96*std(CaConcEnsb,[],2)    CaConcEnsbMean+1.96*std(CaConcEnsb,[],2)];   % 95% Confidence Limits

TimeEnsb = [0:length(CaConcEnsbMean)-1]';       % Define Time (Index) Vector for Ensemble Plots

% PLOT THE DATA, MEAN, AND CONFIDENCE INTERVALS 
figure(7)
hDataLines = plot(TimeEnsb*Ts, CaConcEnsb, '-b');
hold on
hMeanLine = plot(TimeEnsb*Ts, CaConcEnsbMean, '-r', 'LineWidth',2.5);
hCI95Lines = plot(TimeEnsb*Ts, CaCondErr95,   '-r', 'LineWidth',1.5);
hold off
%           (CONSTRUCT ‘LEGEND’ INFORMATION)
hDataLinesGroup = hggroup('DisplayName','Data');
hMeanGroup = hggroup('DisplayName','Mean');
hCI95LinesGroup = hggroup('DisplayName','95% Confidence Limits');

set(hDataLines, 'Parent', hDataLinesGroup)
set(hMeanLine, 'Parent', hMeanGroup)
set(hCI95Lines, 'Parent', hCI95LinesGroup)

set(get(get(hDataLinesGroup, 'Annotation'), 'LegendInformation'), 'IconDisplayStyle', 'on');
set(get(get(hMeanGroup, 'Annotation'), 'LegendInformation'), 'IconDisplayStyle', 'on');
set(get(get(hCI95LinesGroup, 'Annotation'), 'LegendInformation'), 'IconDisplayStyle', 'on');
%           (DISPLAY LEGEND & DEFINE TITLES & AXIS LABELS)
hLgnd7 = legend('show');
set(hLgnd7, 'FontName',Grafont)
set(gca, 'FontName',Grafont, 'FontSize',10)
title('\bfZebrafish Heart Calcium Concentrations\rm', 'FontName',Grafont, 'FontSize',14)
xlabel('Time (s)', 'FontName',Grafont, 'FontSize',12)
ylabel('[Ca^{+2}]  m \itM\rm', 'FontName',Grafont, 'FontSize',12)
axis([0  0.525    -50  550])
orient landscape
grid

figure(8)
contourf(CaConcEnsb)

% ===========================  END CODE  ===========================

Answer 4

[Wayne King]

Michael,

That is helpful. You can see by the following that the dominant oscillation occurs at approximately 2.15 Hz. Let ts be the data vector.

 ts = ts-mean(ts);
 tsdft = fft(ts);
 Fs = 63.2;
 freq = 0:Fs/length(ts):Fs/2;
 plot(freq,abs(tsdft(1:length(ts)/2+1)),'b')
 ylabel('Magnitude'); xlabel('Hz');

The other peaks are harmonics of that dominant 2.15 Hz oscillation. Keep in mind however, that the appearance of harmonics indicates that the oscillation is not perfect (as I would expect). I mean it is a natural oscillatory phenomenon, so a perfect sinusoid is not going to be a perfect match.