Comparing curves to an original signal

I'm trying to compare two curves. They represent the same situation. One curve is the original system and a second this system is undergoing some changes. I would like to compare them. What changes made to the second curve affect the first. And knowing these variations, is it possible for me to go back to the original curve? My curves are experimental data.

I thought about choosing an equation or model that best describes my curve and fitting it using curve fitting or optimization. But I still can't think of something to compare them

[Amanda Botelho Amaral]

these are the curves. Yellow is the original and purple is modified.

[Mathieu NOE]

hello Amanda

may I ask you what is changing indeed ? the system or the input(s) on it ? are the curve a "transfer" characteristic or the output of the system ? what does the x and y axis mean ?

one curve seems so clinical like it' a model and the other one is experimental but it's not what your saying...

now I also don't really know if the purpose is to simply define a metric showing how close the curves are (like a distance or a correlation / fit coefficient) or if it's about trying to find a description how a model can deliver one or the other curve if we know what parameter(s) come into play between the two situations

hope I have made the problem more confusing...:)

[Amanda Botelho Amaral]

Hello Mathieu

The curves are responses to an impedance ground. The yellow curve is ground alone, with nothing connected to it. The purple curve is the response of that same ground when I connect with the tower. So this is the difference, one represents the isolated system and the other the system at work. I would like to get the curve in yellow, starting from the purple one.

just that: "trying to find a description how a model can deliver one or the other curve if we know what parameter(s) come into play between the two situations"

I'm totally lost... :(

[Mathieu NOE]

sorry for my last sentence it was just making the whole think more confusing - forget it

just now I don't really see what you mean by "I would like to get the curve in yellow, starting from the purple one."

do you mean reduce the initial vertical gap so that both curves are overlaid ?

[Amanda Botelho Amaral]

I think I confused you, sorry!

Something similar.... I would like to know a variable, or some element that I can compare.

For example, I have the equation of a curve (Linear), if I change the angular coefficient of this line, I have a new curve. So if I compare the two, can I know what the angular coefficient of each one is? I thought of something like that....

Sorry for my English, I'm learning!

Answer 1

[Jon]

In terms of frequency response of a transfer function it looks like you could make the orange curve move towards the purple curve by adjusting the gain of the transfer function to change the low frequency value. You could provide a pole in the denominator whose value would adjust the high frequency roll off.

[Amanda Botelho Amaral]

Interesting.

Makes a lot of sense. Any tips on how I can do this?

[Jon]

Do you have the control system toolbox?

If so you could for example define a transfer function using its zeros, poles and gains using something like this as as starting point

% define frequencies where curve will break
fz = 1e-7; % zeros, bring magnitude up
fp = [1e-8;1e-5;1.1e-5]; % poles bring magnitude down

% define steady state gain
ssgain = 40

% convert to radians/sec
zeroFit = 2*pi*fz;
poleFit = 2*pi*fp;

% calculate gain
gain = ssgain*prod(poleFit)/prod(zeroFit);

% define transfer function
tf = zpk(zeroFit,poleFit,gain)

% calculate the frequency response 
f = logspace(-9,-4,1000);
H = freqresp(tf,f,'Hz');

% data is returned in a m by n by numpoints array 
% we just want a vector of complex values for SISO system
H = squeeze(H);

% plot the frequency response using semilog plot
semilogx(f,abs(H))

[Amanda Botelho Amaral]

interesting!

I don't know if I have it, but I'll arrange it.

How did you get these poles and zeros?

[Jon]

I will tell you more in a moment but I just saw a glitch in the code I gave you, I named the vector of zero locations zeros, but that is of course the built in MATLAB function which gives zeros. I will edit my answer to correct that.

[Jon]

Each real pole provides a roll off (decreasing amplitude as a function of frequency) starting at a break frequency of approximately the radian locaton of the pole.

Each real zero provide a roll up (increasing amplitude as a function of frequency) starting a a break frequency of approximately the radina location of the zero.

So I placed on pole at 2*pi*1e-8 corresponding to where your function starts to roll off. I also initially place one pole at 2*pi*1e-5 where your function again rolls off. The decrease didn't look rapid enough so I added another pole close by at 2*pi*1.1e-5. To make the flat region, I placed a zero at 2*pi*1e-7 where the flat region begins. The increasing amplitude of this zero cancels the decrease from the low frequency pole, so the net is no increase with frequency in this region, giving the flat spot.

You can adjust those to get a better approximation, I just did it very roughly to show the idea. You can also add more to get steeper slopes.

The gain adjusts the low frequency amplitude

Usually transfer function frequency responses are plotted on log-log scales (Bode Plots) this gives a linear decrease and increase from each real pole, or zero respectively. Might be good to plot it that way.

[Jon]

You can also do this without the Control System Toolbox, just evaluating the complex freqency response directly as the ratio of two polynomial evaluated for a range of purely imaginary frequencies (along the jw axis) Here's code to do that. Some is a repeat of earlier approach using the toolbox, but I thought I'd put it all here so you have it all in one place:

% define frequencies where curve will break
fz = 1e-7; % zeros, bring magnitude up
fp = [1e-8;1e-5;1.1e-5]; % poles bring magnitude down
% define steady state gain
ssgain = 40
% convert to radians/sec
zeroFit = 2*pi*fz;
poleFit = 2*pi*fp;
% calculate gain
gain = ssgain*prod(poleFit)/prod(zeroFit);
% define transfer function
num = poly(zeroFit); % polynomial with roots at zero locaions
den = poly(poleFit); % polynomial with roots at pole locations

% define anonymous function giving the frequency response
fresp = @(s) gain*polyval(num,s)./polyval(den,s);

% calculate the frequency response
w = 2*pi*logspace(-9,-4,1000)
jw = complex(0,w);
H = fresp(jw);

% plot the frequency response magnitude using semilog plot
f = w/2*pi;
semilogx(f,abs(H))

[Jon]

I also will note that if you have the System Identification Toolbox, and maybe the Curve Fitting Toolbox (I'm not familiar with that one) you can fit functions (e.g. transfer functions) to experimental frequency response data. The approach above is simple and may be enough for you though.

[Amanda Botelho Amaral]

Thank you so much for the tip!

helped me a lot

Helped me a lot! I ended up commenting in the comment above about another question. Would you help me?

[Jon]

Hi,

It looks like you were able to solve your second problem.

For the future, it would be good to start a new question if you have an new unrelated, or only distant related issue. That way the answer threads stay clean, and later people can find answers if they have the same question.

Answer 2

[Mathieu NOE]

hello again

so the yellow curve can be shifted to match the purple one by using a linear equation like

y2_shifted = = y2*a+b;

the code provide some initail values for a and b but this could be refined by optimisation

a = ratio2*0.963

b = 1.4

the values 0.963 and 1.4 are my manual fine tuning , but again we could implement a more elegnat and automatic tuning.

see the results now : the yellow has been shifted as close as possible to the purple curve

I just considered that the last portion of the purple curve where the signal sinks would not be taken into account in this procedure...

code :

% extract_data_from_figures
data = extract_data_from_figures('ATE.fig');
x = data.Y(:,1);
y1 = data.Y(:,2);
y2 = data.Y(:,3);

% save amanda.mat x y1 y2 % save in case 

% y ratio between the two curves - measured at first x position
ratio2 = y1(1)./y2(1);
y2_shifted = y2*ratio2*0.963+1.4; 

figure(1),semilogx(x,y1,x,y2,x,y2_shifted,'r');
legend('y1','y2','y2*ratio2*0.963+1.4');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = extract_data_from_figures(filename)
%%
% 
% Input  : File name <filename.fig> 
%          with multiple plots in a single file
% Output : struct data
%        : data.names  contains names of the display object Yvalues
%        : data.Y      contains the actual plot values withthe first column
%                      containing the x-values
%
% Written by Chetanya Puri, 2019
% Last Modified: Nov 6, 2019
%

fig = openfig(filename); % Open figure and assign it to fig object
dataObjs = findobj(fig,'-property','YData'); % Find all graphic objects with YData, in our case line values
xval        = dataObjs(1).XData; % Find the X-axis value
Ymat = [xval(:)]; % Create a matrix with first column of x values
for i=1:length(dataObjs)
    legend_name{i,1} = dataObjs(i).DisplayName;
    yval        = dataObjs(i).YData;
    Ymat = [Ymat yval(:)]; % Keep appending column vectors
end
close(fig); % close the figure
data.names = ['X';legend_name]; 
data.Y = Ymat;
end

[Amanda Botelho Amaral]

Good idea!

I was taking the polynomials for each curve and comparing each of the coefficients to see which one I could standardize for change.

But your way is much simpler!

I will test here this way too!

[Amanda Botelho Amaral]

Thank you so much for the tip! It was a great help

Could you help me with one more thing? I'm having problems with ifft. When I plot the curve, it looks right. However, when zooming in, I see that it is duplicated at the bottom. I don't know what I did wrong.

Here's my code and data.

y=ZZ;

d_vec = (1:(n))./0.2;

t_vec = (1./freq).*d_vec;

time = ifft(y, [], 2);

[Mathieu NOE]

hello

I simply executed this very slightly changed code and I got no such effect

I plot in log x scale so the signal trend is easier in my opinion

load('freq.mat');
load('Z.mat');
n = length(freq);
d_vec = (1:(n))./0.2;
t_vec = (1./freq).*d_vec;
ytime = ifft(ZZ, [], 2);

semilogx(t_vec,ytime);

[Amanda Botelho Amaral]

I found the problem.

I forgot to put the real part of the curve to apply ifft.