curve fitting of a nyquist plot
Show older comments
hi,
I have impedance data points in a nyquist plot:
Anybody an idea how to fit a curve to that data? The data is stored in a array with complex impedance values. I tried a few things, so far without success...
Highly appreciate your support!
2 Comments
Walter Roberson
on 29 Aug 2023
When I follow the plot by eye, I get the impression that for any given Re value that the Im value is either positive or negative, but not both. I do not notice any locations in which the same Re has both a positive and negative Im reading.
Am I mistaken about that, or am I correct but it is an artifact, or is it an actual property of the system?
The question becomes whether there is a continuous curve that is vaguely oval-ish, or if instead it is rather discontinuous and we have to find something in the single parameter Re that jumps back and forth between + and -
Ruben Messner
on 29 Aug 2023
Accepted Answer
Star Strider
on 29 Aug 2023
1 vote
I doubt that there is any direct way to fit it, for example to a nyquist plot. You would have to estimate the system first. The data in the plot are apparently from a Fourier transform, so you have frequency response data. If frequency response data are all you have (rather than the original time-domain data), then use the System Identification Toolbox to identify it, beginning with the idfrd function and going from there, probably ssest since it is the most robust in my experience, although tfest is also an option. Use the compare function to understand graphically how well the estimated system fits the data. If you have the original time-domain data, start with the iddata function.
15 Comments
Ruben Messner
on 29 Aug 2023
Ruben Messner
on 29 Aug 2023
I tried following code:
estimateData = iddata(transpose(voltageFilteredAC(2,:)),transpose(currentFilteredAC(2,:)),Ts(2));
sys = ssest(estimateData, 10);
compare(estimateData, sys);
nyquist(sys);
comparison looks like this:
I used different sys orders, but this is the best. Not good enough I guess.
my manually calculated nyquist:
via nyquist function:
Any idea how to improve the model?
Star Strider
on 29 Aug 2023
First, using time-domain data is preferable, especially if you have the input signal as well as the output. Having the input data would be important in this instance because of the initial transient.
I do not have your data, so I cannot experiment with it. One option would be to increase the system order until you get a decent fit on the compare plot. (I also do not know what your particular system is, and while ssest is usually the best option, experiment with tfest if you cannot get a good fit with ssest and an acceptable order.)
Slightly editing the data to begin wiith the start of the initial square pulse rather than the initial transient may also improve the estimate. (That would have the start time at about 0.5 seconds.) Just send that part of the data (input and output) to iddata first (as a separate iddata structure), estimate that, then include the initial transient later (with the iddata structure you are using here) if you want to.
Ruben Messner
on 29 Aug 2023
Star Strider
on 29 Aug 2023
Thank you!
I need to run this elsewhere (currently MATLAB Online), since it takes more than the allotted 55 seconds to run here. I will keep working on this.
Also, although I assume the data are sampled, the system Identification Toolbox is running this as a continuous-time model. I am not certain if specifying a sampling interval (assumed to be uniform over the entire data set) would be helpful, however it might be worth experimenting with a sampled data set as well. If the data are not uniformly sampled, specify a sampling frequency and then use the Signal Processing Toolbox resample function to do that. The syntax I usually use is:
Fs = 1/round(mean(diff(t)))
[sr,tr] = resample(s,t,Fs);
where ‘t’ is the original time vector, ‘s’ is the original signal vector, and ‘tr’ and ‘sr’ are their uniformly-resampled versions.
Meanwhile, I will keep working with the provided data. (It could take a while.)
@Star Strider I'm thinking that you meant
Fs = round(1/mean(diff(t)))
Instead of
Fs = 1/round(mean(diff(t)))
An example of why I think this
t = 0:0.0015:100;
Fs = 1/round(mean(diff(t)))
Fs = round(1/mean(diff(t)))
Star Strider
on 29 Aug 2023
Yes!
Ruben Messner
on 30 Aug 2023
Star Strider
on 30 Aug 2023
I discovered iddata2timetable (new in R2023a) and was able to determine the sampling interval from the available times. Adding that as an additional name-value pair argument to the estimation functions helped considerably with the accuracy of the estimation, and significantly reduced the time for the estimations to converge (well within the 55-second limit to run them here, although running them in MATLAB Online for convenience). I am still not getting what I consider to be good fits to the data, so I am going to see if additional model options will provide better fits (n4sid, polynomial and ARIMA models, among others).
I worked on this for several hours yesterday. This is not a trivial problem. Since you accepted my Answer, I intend to keep working on it until I have a workable solution.
Ruben Messner
on 31 Aug 2023
I am not happy with these results (I expected better), however they are the best I can come up with. This approach makes the code relatively straightforward to experiment with. The one improvement I added was to add the sampling interval to the estimate functions. This speeds up the code considerably, and provides much better results than considering it as a continuous model.
I included a Fourier transform of the signals, demonstrating band-limited noise. A lowpass filter with a cutoff of 2 Hz should elliminate most of it. Teh best way to do that would be:
uy = lowpass(TT{:,[1 2]}, 2, 1/Ts, 'ImpulseResponse','iir');
u1 = uy(:,1);
y1 = uy(:,2);
You have the relevant signals, so you can create a timetable from the data to replace ‘TT’ here, since the iddata2timetable function was introduced in R2023a. (I needed to do that to see what the signals were, and to get the sampling interval, ‘Ts’.) I printed a segment of it here so you can create one to match it. You can either create it directly, or first create a table and then use the table2timetable function.
These appear to be the only models (I considered several and read the documentation on them) that are appropriate for your data —
LD = load('pemfcDataSet.mat')
estimateData = LD.estimateData;
TT = iddata2timetable(estimateData); % New In R2023
TT(1:5,:)
% TT.Properties
% format long
% t_stats = [mean(diff(TT.Time)); mode(diff(TT.Time)); std(diff(TT.Time)); min(diff(TT.Time)); max(diff(TT.Time))]
% format short
[h,m,s] = hms(TT.Time);
Ts = mean(diff(s))
[FTs1(:,1),Fv] = FFT1(TT{:,1}, s);
[FTs1(:,2),Fv] = FFT1(TT{:,2}, s);
figure
semilogx(Fv, mag2db(abs(FTs1)*2))
grid
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
legend('u_1','y_1', 'Location','best')
xline(2, '--k', '2 Hz')
figure
plot(TT.Time, TT{:,[1 2]})
grid
legend('Input','Output', 'Location','best')
% estimateData = iddata(transpose(voltageFilteredAC(2,:)),transpose(currentFilteredAC(2,:)),Ts(2));
% OPTIONS: 'ss', 'tf', 'armax'
Model = 'ss' % Select Model To Run
switch Model
case 'ss'
tic
sys = ssest(estimateData, 6, 'Ts',Ts)
toc
% % % Order 6 works (sort of)
case 'tf'
tic
sys = tfest(estimateData, 6, 'Ts',Ts)
toc
% % % Order 4 just models the noise
case 'armax'
tic
sys = armax(TT.u1, TT.y1, [3 3 3 1], 'Ts',Ts)
toc
end
figure
compare(estimateData, sys)
figure
nyquist(sys)
function [FTs1,Fv] = FFT1(s,t)
s = s(:);
t = t(:);
L = numel(t);
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv);
end
I will be glad to contnue working on this as long as necessary to get a good result, however this is the best I can do for now. If you have any further ideas, please share them.
.
Ruben Messner
on 1 Sep 2023
Thank you!
Together, no, because the sampling frequency is an important parameter. Models with different sampling frequencies would have different parameters, however the model order can be the same.
An example with a completely synthetic signal demonstrates that the same essential signal (the lengths are different) with different sampling frequencies and the same model order will give models different parameters —
Fs1 = 1000;
L = 500;
t1 = linspace(0, Fs1*L, Fs1*L+1).'/Fs1;
u1 = [0; 1; zeros(L*Fs1-1,1)];
y1 = exp(-0.005*t1) .* sin(2*pi*t1*0.05);
data1 = iddata(y1, u1, 1/Fs1);
figure
plot(t1, y1)
grid
sys1 = ssest(data1,3)
figure
compare(data1, sys1)
Fs2 = 500;
t2 = linspace(0, Fs2*L, Fs2*L+1).'/Fs2;
u2 = [0; 1; zeros(L*Fs2-1,1)];
y2 = exp(-0.005*t2) .* sin(2*pi*t2*0.05);
data2 = iddata(y2, u2, 1/Fs2);
figure
plot(t2, y2)
grid
sys2 = ssest(data2, 3)
figure
compare(data2, sys2)
The pole-zero placements and other characteristics of the two models will be different as werll.
So, with that in mind, you can do something like this with your different inputs and outputs. This uses ssest, however you will need to determine the correct procedure for your data.
.
Ruben Messner
on 6 Sep 2023
thx for that demonstration!
Yeah so merging multiple iddata works only if the sampling frequency is the same I guess...
I also experimented with the system estimation and found out, that removing the offset of the current signal helps a lot. Here is an example:
I calculated the nyquist plot for all dataSets and that's what they look like:
The numbers don't look too bad. Some of them are pretty similar to my other calculations. Maybe a dumb question: Is the nyquist always symmetric when I use nyquist(sys)?
Star Strider
on 6 Sep 2023
As always, my pleasure!
Removing the D-C offset eliminates a significant problem in the estimation. I usually do that by subtracting the mean of the signal. (I did not do that because they were your data and you need to make those decisions.) There is absolutely nothing wrong with doing that providing it does what you want.
I believe the Nyquist diagram should always be symmetric, at least with respect to the real axis. (It is of course not going to be symmetric with respect to the imaginary axis.)
I am happy that it is working and matching your other calculations! I was concerned that it might not work that well, considering that fitting them is the objective here (although that cannot be done directly, at lest in my experience).
More Answers (0)
Categories
Find more on Analyze Data in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
