Phase Noise - Kasdin Algorithm

Hi everyone I'm currently working on phase noise 1/f^alpha. Searching for some code, I didn't found a clear example of implementation. So, I've decided to create my own based on Kasdin. Follow my code:

function [X] = PHkasdin(npts, Qd, alfa)

% npts = number of points % Qd = phase noise variance % alfa = kind of phase noise

    nn = 2*npts;
    ha = alfa/2;
    Qd = sqrt(Qd);
    hfa = zeros(nn,1);
    wfa = zeros(nn,1);
    X = zeros(npts,1);    
    hfa(1) = 1;
    wfa = Qd.*randn(size(wfa));    
    for i = 2: npts     %gera os coeficientes hk
        hfa(i) = hfa(i-1)*( (ha + (i-2))/(i-1) );
%         wfa(i) = Qd*randn();
    end    
    HFA = fft(hfa,nn);
    WFA = fft(wfa,nn);    
    WFA(1) = WFA(1)*HFA(1);
    WFA(2) = WFA(2)*HFA(2);    
    for i = 3: 2: nn
        wr = WFA(i);
        wi = WFA(i+1);
        WFA(i) = wr*HFA(i) - wi*HFA(i+1);
        WFA(i+1) = wr*HFA(i+1) + wi*HFA(i);
    end     
    wfa = (ifft(WFA,npts));    
    for i = 1: npts
        X(i) = X(i) + wfa(i)/npts;
    end
end

In order to get the graph, I wrote the following code:

clear, clc

% signal parameters Fs = 1e6; % sampling frequency, T = 5; % signal duration, s N = round(Fs*T); % number of samples

alfa = 3; % sigma3 = 2.48e-8; sigma3 = 1.26e-4; % just a variance

% noise generation x3 = PHkasdin(N, sigma3 , alfa);

% calculate the PSD winlen = max(size(x3)); na = 8; window = hanning(floor(winlen/na)); % winlen = Fs/2; % window = hanning(winlen, 'periodic');

[Pxx, f] = pwelch((x3), window, 0, [], Fs, 'twosided'); PxxdB = 10*log10(Pxx);

% plot the PSD figure(1) semilogx(f, PxxdB, 'r', 'LineWidth', 2) grid on hold on;

The resulting graph doesn't have a -30dB/decade slope, in fact have a -20dB/decade.

Could someone help me on this?

Thank you

[Mitsuo Sakamoto]

Does the Kasdin's algorithm support alfa > 2?

[Tiago Magalhaes]

Yes, it works for any alfa integer.