Digital signal processing - average fit of periodic signal?

6 次查看（过去 30 天）
Michael
Michael 2012-8-22
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.

请先登录，再回答此问题。

回答（4 个）

Wayne King
Wayne King 2012-8-22
编辑：Wayne King 2012-8-22
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');
  3 个评论
显示 1更早的评论
Wayne King
Wayne King 2012-8-22
编辑：Wayne King 2012-8-22
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
Michael 2012-8-23
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).

请先登录，再进行评论。

Jürgen
Jürgen 2012-8-22
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
  2 个评论
Michael
Michael 2012-8-22
编辑：Michael 2012-8-22
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
Jürgen 2012-8-23
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

请先登录，再进行评论。

Star Strider
Star Strider 2012-8-22
编辑：Star Strider 2012-8-22
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.)
  2 个评论
Michael
Michael 2012-8-22
Thanks for the tip Star Strider, I will give the Peak Fitter a try!
Star Strider
Star Strider 2012-8-26
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 ===========================

请先登录，再进行评论。

Wayne King
Wayne King 2012-8-23
编辑：Wayne King 2012-8-23
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.

类别

在 Help Center 和 File Exchange 中查找有关 Design Condition Indicators Interactively 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by