Моделирование фазового шума во временной области, Как сгенерировать синус с фазовым шумом в MatLab/Octave/SciLab??? Опции

[:-)]

Стало интересно как можно сгенерировать неущее колебание с наложенным фазовым шумом для моделирования система ФАПЧ. Google/yandex неохотно подсказывают как это сделать.

Вот один из найденных вариантов: http://www.dsprelated.com/showcode/246.php

Если я правильно понял идею, то сначала генерируется белый шум, от него берется БПФ, потом на полученный спектр накладывается требуемая форма спектра фазовых шумов и в конце берется ОБПФ.

Неясно как это использовать, если хочется дополнительно наложить меняющийся во времени Доплер на несущую. Есть ли ещё какие способы моделирования фазового шума во временной области?

[Tano]

Покритикуйте программу для расчёта SNR на выходе АЦП из-за джиттера.
Форма сигнала может быть любой.

Код

function Jitter_SNR
    % -- Рачёт зависимости отношения сигнал/шум на выходе АЦП c бесконечной
    % -- разрядностью, от величины джиттера клока, тактирующего АЦП для
    % -- сигнала любой формы.
    % в функции SignaL(t) задаём любой сигнал
    % delta_JIT - шаг джиттера [сек]
    delta_JIT = 1/10^10;                % шаг джиттера [сек]
    imax = 200;                         % количество точек на графике
    % область изменения джиттера: delta_JIT...imax * delta_JIT
    Nmax = 1000;                         % Length of signal
    Fs = 200*10^3;                      % Sampling frequency [Гц]
    T = 1/Fs;                           % Sample time [сек]
    jmax = 100;    
    % считаем напряжение шума, вызванного джиттером с шагом  delta_JIT
    for i = 1:imax                  
        gitter = delta_JIT * i;             % задаём все значения джиттера
        for j = 1:jmax                          
            t = ((j-1)/jmax)*T + (0:Nmax-1)*T;  % Time vector [сек]
            Sref = SignaL(t);               % задаём сигнал на входе АЦП
                                            % он может быть любой
            RMS_Sref(j) = sqrt(sum(Sref.*Sref)/Nmax); % RMS сигнала                                  
            rr = randn(1, Nmax);         % случайный вектор с нормальным
                                         % распределением и единичной дисперсией
            Err_t = gitter*rr;           % вектор джиттера, соответствующий
                                         % этому распределению (в сек)
            tgit = t + Err_t;            % временной вектор с учётом  джиттера
                                         % ( этим тактруем АЦП )
            Sizm = SignaL(tgit);                % сигнал на выходе АЦП
                                                % тактируемым с джиттером
            errSIG = Sref - Sizm;               % вектор ошибки, вызванный джиттером
            U_err(j) = sqrt(sum(errSIG.*errSIG)/Nmax); % RMS вектора ошибки,
                                                       % вызванной джиттером
        end
        RMS(i) = sum(RMS_Sref)/(jmax);
        RMS_error(i) = sum(U_err)/(jmax);    %
    end
    RMS_S = sum(RMS)/imax                  % RMS сигнала   chirp
    DB_error = 20*log10(RMS_error/RMS_S);   % перевод в дб
    x = 1:imax;
    xg = (20/imax)*x;
    figure(1);plot(xg, DB_error); grid;
    title(' SNR [в дб]')
    xlabel('джиттер [нсек]')
end
% --- Синус сигнал -----------------
% function S = SignaL(t)
%     Fsig = 1*10^3;                      % Signal frequency   [Гц] для синуса
%     S = sin(2*pi*Fsig*t);
% end
% --- ЛЧМ сигнал -----------------
function S = SignaL(t)
    Max_t = max(t);
    Fsig = 36*10^3;                      % Signal frequency   [Гц] для синуса
    S = chirp(t, Fsig-Fsig/5, Max_t, Fsig+Fsig/5);     % ЛЧМ сигнал
end

С уважением, Тано