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

Question

Ulrik William Nash 2017-5-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/342235-steps-for-estimating-alpha-in-1-f-alpha-noise-using-ols

回答： Nachiket Katakkar 2017-6-1

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?

3 个评论
显示 1更早的评论隐藏 1更早的评论

Ulrik William Nash 2017-5-28

编辑：John D'Errico 2017-5-28

在 MATLAB Online 中打开

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 2017-5-28

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 2017-6-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/342235-steps-for-estimating-alpha-in-1-f-alpha-noise-using-ols#answer_269242

在 MATLAB Online 中打开

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.