6

I have two vectors that represents a function f(x), and another vector f(ax+b) i.e. a scaled and shifted version of f(x). I would like to find the best scale and shift factors.

*best - by means of least squares error , maximum likelihood, etc.

any ideas?

for example:

f1 = [0;0.450541598502498;0.0838213779969326;0.228976968716819;0.91333736150167;0.152378018969223;0.825816977489547;0.538342435260057;0.996134716626885;0.0781755287531837;0.442678269775446;0];
f2 = [-0.029171964726699;-0.0278570165494982;0.0331454732535324;0.187656956432487;0.358856370923984;0.449974662483267;0.391341738643094;0.244800719791534;0.111797007617227;0.0721767235173722;0.0854437239807415;0.143888234591602;0.251750993723227;0.478953530572365;0.748209818420035;0.908044924557262;0.811960826711455;0.512568916956487;0.22669198638799;0.168136111568694;0.365578085161896;0.644996661336714;0.823562159983554;0.792812945867018;0.656803251999341;0.545799498053254;0.587013303815021;0.777464637372241;0.962722388208354;0.980537136457874;0.734416947254272;0.375435649393553;0.106489547770962;0.0892376361668696;0.242467741982851;0.40610516900965;0.427497319032133;0.301874099075184;0.128396341665384;0.00246347624097456;-0.0322120242872125]

simulated example: scale is 2+7/5, shift 42/55, resample bicubic

*note that f(x) may be irreversible...

Thanks,

Ohad

Mercury
  • 1,886
  • 5
  • 25
  • 44
  • Just to confirm: You don't mean: `a*f(x)+b`? And assuming you don't mean that: You don't know the value of x or what kind of function f is? – Dennis Jaheruddin Nov 26 '12 at 16:59
  • @Rasman: f2(x) = f1(a*x+b). I've given and example for f1 and f2 when a = 3.4 and b = 7/11 – Mercury Nov 27 '12 at 07:27
  • @Mercury Then how can you have more data points in f2 then in f1? You need to have at least some level of consistency – Rasman Nov 27 '12 at 15:22
  • @Rasman You can assume that the total size for the both is the size (i.e. the size of f2), and all the missing elements in f1 are zeros. what I wanted to point is that unlike the "common" problem of f2(i)=a*f1(i)+b, in this case f1(i) and f2(i) are not necessarily dependent in any way – Mercury Nov 28 '12 at 06:58
  • It appears that more information is required, do you really just have 2 vectors and want to use them to discover how they were created? Because in that case it is impossible to find a best solution. – Dennis Jaheruddin Jan 21 '13 at 10:20
  • Okay, there is one thing I still don't quite get. f1(x) = f(x); f2(x)=f(a*x+b) this is fine. But the two vectors you got, the x they would be created with are not necessarily of the same length (or similar at all)? If we could assume the values were created by a common vector x i guess this would be possible, but if we are dealing with 2 different vectors x that are independent, I don't think there is a solution. – Dani Gehtdichnixan Jan 24 '13 at 15:25

5 Answers5

6

For each f(x), take the absolute value of f(x) and normalize it such that it can be considered a probability mass function over its support. Calculate the expected value E[x] and variance of Var[x]. Then, we have that

  E[a x + b] = a E[x] + b
Var[a x + b] = a^2 Var[x]

Use the above equations and the known values of E[x] and Var[x] to calculate a and b. Taking your values of f1 and f2 from your example, the following Octave script performs this procedure:

% Octave script
% f1, f2 are defined as given in your example
f1 = [zeros(length(f2) - length(f1), 1); f1];

save_f1 = f1; save_f2 = f2;

f1 = abs( f1 ); f2 = abs( f2 );
f1 = f1 ./ sum( f1 ); f2 = f2 ./ sum( f2 );

mean = @(x)sum(((1:length(x))' .* x));
var = @(x)sum((((1:length(x))'-mean(x)).^2) .* x);

m1 = mean(f1); m2 = mean(f2);
v1 = var(f1); v2 = var(f2)

a = sqrt( v2 / v1 ); b = m2 - a * m1;

plot( a .* (1:length( save_f1 )) + b, save_f1, ...
      1:length( save_f2 ), save_f2 );
axis([0 length( save_f1 )];

And the output is enter image description here

jstarr
  • 480
  • 3
  • 4
  • oooh I like this one ... :) – Rody Oldenhuis Jan 27 '13 at 21:24
  • If you replace the absolute value step with squaring, then one could think of finding `a` as answering the question "what (time) compression factor do I need to make the new signal energy equal to the old signal energy?" It probably would not work so well if the zero-padding was replaced with noise. :) – jstarr Jan 28 '13 at 23:10
5

Here's a simple, effective, but perhaps somewhat naive approach.

First make sure you make a generic interpolator through both functions. That way you can evaluate both functions in between the given data points. I used a cubic-splines interpolator, since that seems general enough for the type of smooth functions you provided (and does not require additional toolboxes).

Then you evaluate the source function ("original") at a large number of points. Use this number also as a parameter in an inline function, that takes as input X, where 

X = [a b]

(as in ax+b). For any input X, this inline function will compute

  1. the function values of the target function at the same x-locations, but then scaled and offset by a and b, respectively.

  2. The sum of the squared-differences between the resulting function values, and the ones of the source function you computed earlier.

Use this inline function in fminsearch with some initial estimate (one that you have obtained visually or by via automatic means). For the example you provided, I used a few random ones, which all converged to near-optimal fits.

All of the above in code: 

function s = findScaleOffset

    %% initialize

    f2 = [0;0.450541598502498;0.0838213779969326;0.228976968716819;0.91333736150167;0.152378018969223;0.825816977489547;0.538342435260057;0.996134716626885;0.0781755287531837;0.442678269775446;0];
    f1 = [-0.029171964726699;-0.0278570165494982;0.0331454732535324;0.187656956432487;0.358856370923984;0.449974662483267;0.391341738643094;0.244800719791534;0.111797007617227;0.0721767235173722;0.0854437239807415;0.143888234591602;0.251750993723227;0.478953530572365;0.748209818420035;0.908044924557262;0.811960826711455;0.512568916956487;0.22669198638799;0.168136111568694;0.365578085161896;0.644996661336714;0.823562159983554;0.792812945867018;0.656803251999341;0.545799498053254;0.587013303815021;0.777464637372241;0.962722388208354;0.980537136457874;0.734416947254272;0.375435649393553;0.106489547770962;0.0892376361668696;0.242467741982851;0.40610516900965;0.427497319032133;0.301874099075184;0.128396341665384;0.00246347624097456;-0.0322120242872125];

    figure(1), clf, hold on

    h(1) = subplot(2,1,1); hold on
    plot(f1);
    legend('Original')

    h(2) = subplot(2,1,2); hold on
    plot(f2);

    linkaxes(h)
    axis([0 max(length(f1),length(f2)), min(min(f1),min(f2)),max(max(f1),max(f2))])


    %% make cubic interpolators and test points

    pp1 = spline(1:numel(f1), f1);
    pp2 = spline(1:numel(f2), f2);

    maxX = max(numel(f1), numel(f2));
    N  = 100 * maxX;

    x2 = linspace(1, maxX, N);
    y1 = ppval(pp1, x2);

    %% search for parameters

    s = fminsearch(@(X) sum( (y1 - ppval(pp2,X(1)*x2+X(2))).^2 ), [0 0])

    %% plot results

    y2 = ppval( pp2, s(1)*x2+s(2));

    figure(1), hold on
    subplot(2,1,2), hold on    
    plot(x2,y2, 'r')
    legend('before', 'after')



end

Results: 

s =
2.886234493867320e-001    3.734482822175923e-001

the result

Note that this computes the opposite transformation from the one you generated the data with. Reversing the numbers: 

>> 1/s(1) 
ans =    
    3.464721948700991e+000   % seems pretty decent 
>> -s(2)
ans = 
    -3.734482822175923e-001  % hmmm...rather different from 7/11!

(I'm not sure about the 7/11 value you provided; using the exact values you gave to make a plot results in a less accurate approximation to the source function...Are you sure about the 7/11?)

Accuracy can be improved by either

  1. using a different optimizer (fmincon, fminunc, etc.)
  2. demanding a higher accuracy from fminsearch through optimset
  3. having more sample points in both f1 and f2 to improve the quality of the interpolations
  4. Using a better initial estimate

Anyway, this approach is pretty general and gives nice results. It also requires no toolboxes.

It has one major drawback though -- the solution found may not be the global optimizer, e.g., the quality of the outcomes of this method could be quite sensitive to the initial estimate you provide. So, always make a (difference) plot to make sure the final solution is accurate, or if you have a large number of such things to do, compute some sort of quality factor upon which you decide to re-start the optimization with a different initial estimate.

It is of course very possible to use the results of the Fourier+Mellin transforms (as suggested by chaohuang below) as an initial estimate to this method. That might be overkill for the simple example you provide, but I can easily imagine situations where this could indeed be very useful.

Rody Oldenhuis
  • 37,726
  • 7
  • 50
  • 96
  • First of all, you were right about the 7/11, I can't remember what value I used... Regarding your solution: As you said although it is very brute force, it works. There only problem with this solution is that it is time consuming (iterative ,extensive bicubic interpolation). But since computation limits were not stated in the question, this Answer is pretty good :) – Mercury Jan 21 '13 at 21:36
  • @Mercury: Fourier+Mellin are also not cheap :) In any case, I suppose you could only use the `x`-locations of the target function, evaluated in the source function, to compute the squared-errors with. This would indeed reduce computation time and should (as far as I can see) have little impact on the accuracy of the outcome. There is a slight risk of sub-Nyquist ambiguities, but oh well. – Rody Oldenhuis Jan 23 '13 at 08:31
  • @Mercury: By the way, I assumed you have no knowledge of `f1`'s analytical representation. Do you know the general analytical form of `f1`? If so, perhaps we can conjure up a more specialized version for such cases. – Rody Oldenhuis Jan 23 '13 at 08:35
  • @Mercury: Also by the way, I assumed you also don't have knowledge about the specific locations and number of sample points in `f1` and `f2`. Should we assume they are always uncorrelated? If they are correlated in some known way, the problem is of course pretty trivial. – Rody Oldenhuis Jan 23 '13 at 08:39
  • f1 and f2 are samples of a arbitrary function f in a different sampling frequency and offset. f can not be analytically written. Since the scale factor is unknown, they are not correlated, and even if the scale was known, correlation is good up to about half a sample accuracy. – Mercury Jan 23 '13 at 12:17
  • @Mercury: Yeah I thought as much. In that case, stick with the method above :) – Rody Oldenhuis Jan 23 '13 at 12:45
3

For the scale factor a, you can estimate it by computing the ratio of the amplitude spectra of the two signals since the Fourier transform is invariant to shift.

Similarly, you can estimate the shift factor b by using the Mellin transform, which is scale invariant.

chaohuang
  • 3,965
  • 4
  • 27
  • 35
  • 1
    The Fourier transform trick was already tried, but the results are not accurate enough, not even on simulated data, nevertheless on real-life problems. I'm was not very familiar with the Mellin transform, but now I am :). The problem here is the same as in the Fourier transform - inaccuracy. Note that If I would have guessed the scale successively, simple correlation would give me the shift. – Mercury Dec 11 '12 at 10:08
1

Here's a super simple approach to estimate the scale a that works on your example data:

a = length(f2) / length(f1)

This gives 3.4167 which is close to your stated value of 3.4. If that estimate is good enough, you can use correlation to estimate the shift.

I realize that this is not exactly what you asked, but it may be an acceptable alternative depending on the data.

shoelzer
  • 10,648
  • 2
  • 28
  • 49
  • Not true if `f2` data comes from another sampling than `f1`, i.e., `f1` and `f2` have a non-related number of sample points. – Rody Oldenhuis Jan 23 '13 at 08:32
  • @RodyOldenhuis You're right, but that was not stated in the question. – shoelzer Jan 23 '13 at 13:51
  • Yes it was, implicitly; through the example he provides. – Rody Oldenhuis Jan 23 '13 at 14:23
  • 1
    @RodyOldenhuis I disagree. In the example, `f1` and `f2` are sampled over the same interval but at different rates. This makes the number of samples related and that's why `a = length(f2) / length(f1)` is close to the right answer. – shoelzer Jan 23 '13 at 14:47
  • Yes, you're right, sorry. Nevertheless, it now turns out the OP indeed *intended* the data to come from independent sources, with independent sample rates. +1 though for being a smartass :) – Rody Oldenhuis Jan 23 '13 at 15:04
1

Both Rody Oldenhuis and jstarr's answers are correct. I'm adding my own answer just to sum things up, and connect between them. I've messed up Rody's code a little bit and ended up with the following:

function findScaleShift
load f1f2

x0 = [length(f1)/length(f2) 0]; %initial guess, can do better
n=length(f1);
costFunc = @(z) sum((eval_f1(z,f2,n)-f1).^2);
opt.TolFun = eps; 
xopt=fminsearch(costFunc,x0,opt);
f1r=eval_f1(xopt,f2,n);
subplot(211);
plot(1:n,f1,1:n,f1r,'--','linewidth',5)
title(xopt);
subplot(212);
plot(1:n,(f1-f1r).^2);
title('squared error')
end


function y = eval_f1(x,f2,n)
t = maketform('affine',[x(1) 0 x(2); 0 1 0 ; 0 0 1]');
y=imtransform(f2',t,'cubic','xdata',[1 n ],'ydata',[1 1])';
end

This gives zero results: enter image description here This method is accurate but exhaustive and may take some time. Another disadvantage is that it finds only a local minima, and may give false results if initial guess (x0) is far.

On the other hand, jstarr method gave the following results:

xopt = [ 3.49655562549115         -0.676062367063033]

which is 10% deviation from the correct answer. Pretty fast solution, but not as accurate as I requested, but still should be noted. I think in order to get the best results jstarr method should be used as an initial guess for the method purposed by Rody, giving an accurate solution.

Ohad

Mercury
  • 1,886
  • 5
  • 25
  • 44