Steps for estimating alpha in 1/f^alpha noise using OLS

I have a signal (the vector "pink"), which I know is pink noise, 1/f^alfa, where alpha = 1. But let's assume alpha is unknown, and I need to estimate alpha. In other words, Let's assume I just have the signal.

I use the following procedure to get into a position to estimate alpha using OLS:

F = abs(fft(pink)); x = log(1:(length(pink)/2)); y = log(F(1:(length(pink)/2)));

Now I run OLS regression:

[r,m,b] = regression(x,y)

where I get

r = -0.5948, m = -0.4915, b = 6.8223.

I would have thought my estimate of alpha was m, but m is not equal to 1. So, what am I doing wrong? Is there an additional step I am forgetting, or is my procedure plain wrong?

[John D'Errico]

I don't know of any function called regression, nor have I ever heard of one. A search of the site does not show any tool by that name, in case there is one in some toolbox. I also do not see anything of that form in the file exchange. My searches might have missed something though. I did look at each of the 5 files named regression.m that I found on the file exchange. None of them seemed to fit that template.

So if you got some code by that name, then ask the author. If it is in some toolbox that I did not find, then tell where it came from. If you wrote it, then how can we know what it does and what the return arguments mean?

[Ulrik William Nash]

Hi John,

I didn't expect that twist.

Here is the code for the particular function:

function [r,m,b] = regression(targets,outputs,flag)
%REGRESSION Linear regression.
%  
%  <a href="matlab:doc regression">regression</a> calculates the linear regression between each element
%   of the network response and the corresponding target.
%  
%  [R,M,B] = <a href="matlab:doc regression">regression</a>(T,Y) takes cell array or matrix targets T and
%  output Y, each with total matrix rows of N, and returns the linear
%  regression for each of the N rows: the regression values R, slopes M,
%  and y-intercepts B.
%
%  <a href="matlab:doc regression">regression</a>(T,Y,'one') returns scalar R, M and B values across all
%  rows of targets and outputs.
%
%  Here a feedforward network is trained and regression performed on its
%  targets and outputs.
%
%    [x,t] = <a href="matlab:doc simplefit_dataset">simplefit_dataset</a>;
%    net = <a href="matlab:doc feedforwardnet">feedforwardnet</a>(10);
%    net = <a href="matlab:doc train">train</a>(net,x,t);
%    y = net(x);
%    [r,m,b] = <a href="matlab:doc regression">regression</a>(t,y)
%
%  See also PLOTREGRESSION

% Copyright 2010 The MathWorks, Inc.

if nargin < 2, error(message('nnet:Args:NotEnough')); end

if iscell(targets), targets = cell2mat(targets); end
if iscell(outputs), outputs = cell2mat(outputs); end

if all(size(targets) ~= size(outputs))
  error(message('nnet:NNData:TYMismatch'))
end

if (nargin >= 3) && ischar(flag) && strcmp(flag,'one')
  targets = targets(:)';
  outputs = outputs(:)';
end

[N,Q] = size(targets);
m = zeros(N,1);
b = zeros(N,1);
r = zeros(N,1);
for i=1:N
  t = targets(i,:);
  y = outputs(i,:);
  ignore = find(isnan(t) | isnan(y));
  t(ignore) = [];
  y(ignore) = [];
  Quse = Q - length(ignore);
  h = [t' ones(size(t'))];
  yt = y';
  rankStatus = warning('off','MATLAB:rankDeficientMatrix');
  rankRestore = onCleanup(@() warning(rankStatus));
  theta = h\yt;
  m(i) = theta(1);
  b(i) = theta(2);
  yn = y - mean(y);
  tn = t - mean(t);
  sty = std(yn);
  stt = std(tn);
  r(i) = yn*tn'/(Quse - 1);
  if (sty~=0)&&(stt~=0)
    r(i) = r(i)/(sty*stt);
  end
end

[Ulrik William Nash]

I just double-checked that the "regression" and "regress" functions compute the same output, which they do, except the first also calculates r.

Answer 1

[Nachiket Katakkar]

The regression function is part of the Neural Network Toolbox:

http://www.mathworks.com/help/nnet/ref/regression.html

The coefficient "m" is the slope of the curve between the inputs "x" and "y". The example on the documentation page shows a value of m equal to 1, however, this would not necessarily be the case for any regression problem.

Consider this example with random data:

>> pink = rand(1,100);
F = abs(fft(pink));
x = log(1:(length(pink)/2));
y = log(F(1:(length(pink)/2)));
[r,m,b] = regression(x,y)
r =

   -0.3448

 m =

   -0.3469

 b =

    1.7326

You can visualize the relationship between "x" and 'y" using:

 >> plotregression(x,y)

Also notice, that in your case the value of "b", the offset of regression is also quite high. The regression function seems to be working as expected, so I would recommend confirming that the algorithm you are using is correct.