function x = one_on_f_noise2(N, alpha, seed)
%
% one_on_f_noise.m, (c) M Roughan, 2018
%
% created: 	Thu Mar 22 2018 
% author:  	M Roughan 
% email:   	matthew.roughan@adelaide.edu.au
% 
% ONE_ON_F_NOISE 
%         generate 1/f^alpha noise using a Fourier technique
%
% INPUTS:
%         N = number of frequencies or equivalently number of time points (in x and y directions)
%         alpha = exponent in the filter
%         
% OUTPUTS:        
%        x = NxN signal 
%         
%
%

if nargin < 2
  alpha = 1.0
end

if nargin == 3
  rng('default');
  rng(seed);
end

% generate the Fourier domain signal
magnitude = randn(N/2, N/2); % Normal (Gaussian) random numbers
phase = 2*pi*rand(N/2, N/2); % uniformly random phases
f_quarter = magnitude .* exp( i*phase ); % complex values for frequency components

% now filter it to get a 1/f spectrum
x = (1:N/2).^alpha;
[X,Y] = meshgrid(x);
% f_quarter = f_quarter ./ max(X, Y);
f_quarter = f_quarter ./ sqrt(X + Y);

% make sure the signal is Hermitian
f_full = [f_quarter,               fliplr(conj(f_quarter));
	  flipud(conj(f_quarter)), flipud(fliplr(f_quarter));
	 ];       

% use to Inverse Fourier transform to go back to the time domain
%    making sure to get a real signal, though that shouldn't be needed
x = real(ifft2(f_full));




