Digital signal processing - average fit of periodic signal?
12 ビュー (過去 30 日間)
古いコメントを表示
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:
![](https://www.mathworks.com/matlabcentral/images/broken_image.png)
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.
0 件のコメント
回答 (4 件)
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 件のコメント
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.
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 件のコメント
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
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 件のコメント
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
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.
0 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!