Plot peaks for each column in a csv then fit a curve.
Ältere Kommentare anzeigen
I am trying to process some data obtained from an oscilloscope. The columns in the csv represent the amplitude of the signal. I expected to see what is in the picture. I want to find the peaks for each signal and then plot the peaks as a function of tau (ms) ranging from 0.1 and up. I'd also like to fit a curve to each line following the exponential function y = a*e(-x/T1) +c, where a is a constant (-2) and T1 is found from the curve. Finding T1 is the ultimate goal of plotting the signals. I have tried a few different methods: here is my most recent attempt and it did not give me anywhere near what I was hoping for:
filename = 'D:\Grad Lab\NMR\Data\T1 Data\compiledT1nolabel.csv';
data = readtable(filename);
y1= table2array(data(:,1));
y2= table2array(data(:,2));
y3= table2array(data(:,3));
y4= table2array(data(:,4));
y5= table2array(data(:,5));
y6= table2array(data(:,6));
y7= table2array(data(:,7));
Fs = 2e-10;
T = 1/Fs;
L = 2500;
t = (0:L-1)*T;
Fn = Fs/2;
Fy1 = fft(y1)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:length(Fv);
figure(1)
plot(Fv, abs(Fy1(Iv))*2)
grid
title("Fourier Transfrom of Mineral Oil Original Signal")
xlabel('Frequency (Hz)')
ylabel('Amplitude')
Fy1dcoc = fft(y1-mean(y1))/L;
figure(2)
plot(Fv, abs(Fy1dcoc(Iv))*2)
grid
title('Fourier Transform of D-C Offset Corrected Signal')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
[Fy1n_max,Iv_max] = max(abs(Fy1dcoc(Iv))*2);
Wp = 2*Fv(Iv_max)/Fn;
Ws = Wp*2;
Rp = 10
Rs = 30;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[b,a] = butter(n,Wn);
[SOS,G] = tf2sos(b,a);
S = filtfilt(SOS,G,y1);
figure(4)
plot(t, y1)
hold on
plot(t, S, '-r', 'LineWidth',1.5)
hold off
grid
legend('Mineral oil', "'S'", 'Location', 'N')
title('Original Signal of Mineral Oil and Uncorrupted Signal (S)')
xlabel('Tau (ms)')
ylabel('Amplitude')
I appreciate any help with this!
Akzeptierte Antwort
There are no visible peaks in the time-domain signal.
Doing the Foureir transform and finding the peaks is easy enough, however I have no idea how to do the plot in the image because I have no idea how to calculate whatever the variables are that it depicts. The independent variable appears to be time.
If your intent is to do a time-frequency analysis, it would be best to use the pspectrum function with the 'spectrogram' option, rather than fft.
data = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1178628/compiledT1.csv', 'VariableNamingRule','preserve')
VN = data.Properties.VariableNames;
Fs = 2e-10;
T = 1/Fs;
L = 2500;
t = (0:L-1)*T;
Fn = Fs/2;
y = table2array(data)
figure
plot(t, y)
grid
xlabel('Time')
ylabel('Amplitude')
legend(VN, 'Location','best')
ydt1 = detrend(y(:,1),1);
figure
plot(t, ydt1)
grid
xlabel('Time')
ylabel('Amplitude')
title(VN{1})
[p,f,t] = pspectrum(ydt1,Fs,'spectrogram');
figure
waterfall(f,t,p')
xlabel('Frequency (Hz)')
ylabel('Time (seconds)')
wtf = gca;
wtf.XDir = 'reverse';
colormap(turbo)
view([30 45])
NFFT = 2^nextpow2(L);
ym = y - mean(y);
FTy = fft(ym, NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTy(Iv,:))*2)
grid
xlabel('Frequency')
ylabel('Magnitude')
legend(VN, 'Location','best')
xlim([0 5E-12])
[pks1,locs1] = findpeaks(abs(FTy(Iv,1))*2);
figure
plot(Fv, abs(FTy(Iv,1))*2, 'DisplayName','Data')
hold on
plot(Fv(locs1), pks1, '^r', 'DisplayName','Peaks')
hold off
grid
xlabel('Frequency')
ylabel('Magnitude')
title(VN{1})
legend()
% legend(VN{1}, 'Location','best')
xlim([0 5E-12])
.
9 Kommentare
Alex
am 3 Nov. 2022
data = readtable(amplitudeT1.xlsx)
VN = data.Properties.VariableNames;
Fs = 1e-10;
T = 1/Fs;
L = 2500;
t = (0:L-1)*T;
Fn = Fs/2;
y = table2array(data);
figure(1)
plot(t, y)
grid
xlabel('Time (ms)')
ylabel('Amplitude')
legend(VN, 'Location', 'best')
ydt1 = detrend(y(:,1),1);
figure(2)
plot(t, ydt1)
grid
xlabel('Time')
ylabel('Amplitude')
title(VN{1})
[p,f,t] = pspectrum(ydt1,Fs,'spectrogram');
figure(3)
waterfall(f,t,p')
xlabel('Frequency (Hz)')
ylabel('Time')
wtf = gca;
wtf.xDir = 'reverse';
colormap(turbo)
view([30 45])
NFFT = 2^nextpow2(L);
ym = y - mean(y);
FTy = fft(ym, NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure(4)
plot(Fv, abs(FTy(Iv,:))*2)
grid
xlabel('Frequency')
ylabel('Magnitude')
legend(VN, 'Location', 'best')
xlim([0 5e-12])
[pks1,locs1] = findpeaks(abs(FTy(Iv,1))*2);
figure(5)
plot(Fv, abs(FTy(Iv,1))*2), 'DisplayName','Data')
hold on
plot(Fv(locs1), pks1, '^r', 'DisplayName', 'Peaks')
hold off
grid
xlabel('Frequency')
ylabel('Magnitude')
title(VN{1})
legend(VN{1}, 'Location', 'best')
xlim([0 5E-12])
Running this it looks great (all the plots), I get the same plots you saw. I attached the incorrect file so I attached the correct one. I think the best thing for me to do would be to just plot each one separately though seeing how busy these plots. For the plot attached I want to grab that highest peak and then every other peak after that. Then if I can fit a curve to those peaks I can find T1 (the relaxation time of the atom.) Thank you!
It took a while to get this to work. This is the best I can do with these data.
Part of the problem is that the magnitude of the ‘t’ vector puts many of the calculations near the limit of floating-point precision for addition and subtraction (the eps value). That is also the reason that it’s not possible to design a functioning digital filter for these data. (I gave that a shot as well.)
Also, the ‘a’ value of -2 will absolutely not work unless the data themsellves are also negative. I let that float here.
Editing these data further, for example starting at the maximum peak, (the fourth peak for the Mineral Oil data) may produce a better fit. I don’t know what these data are or what they represent, so I didn’t edit them. I also didn’t detrend them, although they have a wandering baseline.
The Fourier transform data are interesting. All the data have peaks at either 
or 
, and some have both.
or , and some have both. Try this —
data = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1178803/amplitudeT1.xlsx', 'VariableNamingRule','preserve')
VN = data.Properties.VariableNames;
Fs = 2e-10;
T = 1/Fs;
L = 2500;
t = (0:L-1)*T;
Fn = Fs/2;
y = table2array(data);
figure
plot(t, y)
grid
xlabel('Time')
ylabel('Amplitude')
legend(VN, 'Location','best')
% Lminv = islocalmin(y(:,1), 'MinProminence',1);
% p = polyfit(t(Lminv).', y(Lminv,1), 9);
% ydt1 = y(:,1) - polyval(p,t(:));
% ydt1 = detrend(y(:,1),15);
[pkst1, locst1] = findpeaks(y(:,1), 'MinPeakProminence',4.5);
figure
% plot(t, ydt1)
plot(t, y(:,1))
hold on
plot(t(locst1), pkst1, '^r')
% plot(t(Lminv), y(Lminv,1), '+r')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
title(VN{1})
format longE
p = polyfit(t(locst1), log(pkst1), 1) % Estimate Initial Parameters For Exponential Fit
ypv = exp(polyval(p, t(locst1)));
figure
plot(t(locst1), pkst1,'.b', t(locst1),ypv, '-r')
grid
xlabel('Time')
ylabel('Amplitude')
title('Initial Parameter Estimate Result')
% y = a*e(-x/T1) +c
objfcn = @(b,t) b(1).*exp(-t./b(2)) + b(3);
mdl = fitnlm(t(locst1), pkst1, objfcn, [p(2); -1/(p(1)); 1])
tv = linspace(min(t(locst1)), max(t(locst1)));
[yv,yci] = predict(mdl, tv(:));
figure
hp{1} = plot(t(locst1), pkst1,'.b', 'DisplayName','Data');
hold on
hp{2} = plot(tv,yv,'-r', 'DisplayName','Exponential Fit');
hp{3} = plot(tv,yci,'--r', 'DisplayName','±95% CI');
hold off
grid
xlabel('Time')
ylabel('Amplitude')
title(VN{1})
legend([hp{1};hp{2};hp{3}(1)],'Location','best')
[p,f,t] = pspectrum(y(:,1),Fs,'spectrogram');
figure
waterfall(f,t,p')
xlabel('Frequency (Hz)')
ylabel('Time (seconds)')
wtf = gca;
wtf.XDir = 'reverse';
colormap(turbo)
xlim([0 3E-11])
view([135 45])
NFFT = 2^nextpow2(L);
ym = y - mean(y);
FTy = fft(ym, NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTy(Iv,:))*2)
grid
xlabel('Frequency')
ylabel('Magnitude')
legend(VN, 'Location','best')
xlim([0 5E-12])
[pksf1,locsf1] = findpeaks(abs(FTy(Iv,1))*2);
figure
plot(Fv, abs(FTy(Iv,1))*2, 'DisplayName','Data')
hold on
plot(Fv(locsf1), pksf1, '^r', 'DisplayName','Peaks')
hold off
grid
xlabel('Frequency')
ylabel('Magnitude')
title(VN{1})
legend()
% legend(VN{1}, 'Location','best')
xlim([0 5E-12])
.
Alex
am 3 Nov. 2022
Star Strider
am 3 Nov. 2022
My pleasure!
Alex
am 3 Nov. 2022
Alex
am 4 Nov. 2022
Star Strider
am 4 Nov. 2022
Bearbeitet: Star Strider
am 4 Nov. 2022
If the peaks are supposed to be at 1 ms intervals, and you want to fit from the highest peak to the peaks greater than 2, this works —
data = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1178803/amplitudeT1.xlsx', 'VariableNamingRule','preserve')
VN = data.Properties.VariableNames;
Fs = 2e-10;
T = 1/Fs;
L = 2500;
t = (0:L-1)*T;
Fn = Fs/2;
y = table2array(data);
figure
plot(t, y)
grid
xlabel('Time')
ylabel('Amplitude')
legend(VN, 'Location','best')
[pkst1, locst1] = findpeaks(y(:,1)); % Detect All Peaks
figure
plot(t, y(:,1))
hold on
plot(t(locst1), pkst1, '^r')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
title(VN{1})
[mxpk,mxpkidx] = max(pkst1)
Lvpeaks = (pkst1 <= max(pkst1)) & (pkst1 >= 2.45) & (locst1 >= locst1(mxpkidx)); % Select Peaks
locst1e = locst1(Lvpeaks) % Select ‘locs'
pkst1e = pkst1(Lvpeaks); % Select Corresponding 'pks'
tauv = (1:numel(pkst1e))*1E-3; % Create 'tau' Time Vector
format longE
% p = polyfit(t(locst1), log(pkst1), 1) % Estimate Initial Parameters For Exponential Fit
p = polyfit(tauv, log(pkst1e), 1) % Estimate Initial Parameters For Exponential Fit
ypv = exp(polyval(p, tauv));
figure
plot(tauv, pkst1e,'.b', tauv,ypv, '-r')
grid
xlabel('\tau')
ylabel('Amplitude')
title('Initial Parameter Estimate Result')
% y = a*e(-x/T1) +c
objfcn = @(b,t) b(1).*exp(-t./b(2)) + b(3);
mdl = fitnlm(tauv, pkst1e, objfcn, [p(2); -1/(p(1)); 1])
tv = linspace(min(tauv), max(tauv));
[yv,yci] = predict(mdl, tv(:));
figure
hp{1} = plot(tauv, pkst1e,'.b', 'DisplayName','Data');
hold on
hp{2} = plot(tv,yv,'-r', 'DisplayName','Exponential Fit');
hp{3} = plot(tv,yci,'--r', 'DisplayName','±95% CI');
hold off
grid
xlabel('\tau')
ylabel('Amplitude')
title(VN{1})
legend([hp{1};hp{2};hp{3}(1)],'Location','best')
This is the best I can do. I am not certain how to put ‘Peak 2’ specifically into the set of peaks to be regressed, or eliminate the others.
.
Alex
am 4 Nov. 2022
Star Strider
am 4 Nov. 2022
My pleasure!
I also don’t understand the time values on the x-axis in the figure in your last Comment. In my code, I assigned each peak x-value a serial 0.1 millisecond value. That may not be correct, however that’s how I interpreted ‘So for each peak that is a one tau (0.1ms.)’ although it doesn’t make sense to me.
Weitere Antworten (2)
Image Analyst
am 3 Nov. 2022
0 Stimmen
See attached demo. It will be easy for you to adapt it to your data. Just replace demo data with your actual data.
6 Kommentare
Alex
am 3 Nov. 2022
Image Analyst
am 3 Nov. 2022
The curve does not look like an inverted exponential:
Why do you want to fit the red line to an inverted exponential instead of some sigmoid function, like the rate equation (demo atached)?
Alex
am 3 Nov. 2022
Image Analyst
am 3 Nov. 2022
I just plotted s vs t, and y1 vs t like you did. I'm not sure why the plot I posted above looks so drastically different than @Star Strider's.
Well if s and y1 are your measured values (I doubt it since they're so smooth) and if you say that the theory behind the situation says that it's supposed to be an inverted exponential decay, then use an exponential decay.
I thought that my demo was pretty well commented. Did you have trouble replacing y with s, and x with your t? If you can't do it, I can do it for you. But essentially it's just passing x and y and the theoretical model into fitnlm() - just a few lines of code if you ignore the comments and fancy plotting stuff.
Alex
am 4 Nov. 2022
Image Analyst
am 6 Nov. 2022
There are 7 columns in that workbook. What do the different columns represent? Which column is the x value and which are/is the y value(s)? What column are you trying to fit with an exponential?
Alex
am 4 Nov. 2022
0 Stimmen
Sorry for the confusion (if it's any consolation by going through this I am understanding the experiment more.) The signal in the 'attached figure' (not the picture) that is a signal from an oscilloscope. I drew arrows indicating 'peak 1', 'peak 2' and 'peak 3' (the first couple peaks before the labeled peak 1 are noise and can be deleted.) What is done we subject a magnetic field by emitting a pulse (peak 1) then we flip that field with a second pulse at a time (tau = 0.2ms) (peak 2) then there is an echo (peak 3) at 2 tau = 0.4 (the next echo would be at 0.8ms and so on.) I am after the decaying echo peaks that I drew arrows to the majority of them (they follow a decreasing trend. This trend will allow us to find T1 our relaxation time of when the field reverts back to its normal axis.
the equations have been different in each different paper i have read, however, y = A-B *exp(-tau/T1) is used the most. in the indicates the decaying of the echo as tau increases. The field should level out so I only need to really draw it out to tau 100 ms (or shorter if it levels off before then.)
I hope this helps clarify things a bit!
3 Kommentare
Viewing the figure —
F = openfig(websave('peaksmineraloil','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1180283/peaksmineraloil.fig'));
I’ll edit my latest Comment to use only the tallest peak and then all peaks greater than 2 ...
.
Alex
am 6 Nov. 2022
Star Strider
am 6 Nov. 2022
As always, my pleasure!
Kategorien
Mehr zu Mathematics finden Sie in Hilfe-Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
