模拟随机过程的三角级数法中的一个问题

sy0116

三角级数法可根据给定的自功率谱密度函数模拟各种各态历经的高斯过程，其在星谷胜的《随机振动分析》中有详细介绍，其中一种实现形式可用如下公式给出：

其中各谐波项的频率fn由下式给出：
               （1）
其中：
                                                          （2）
按此计算，谐波项的最低频率为 ，最高为
我的问题是：上面（1）式中为何要用n减去1/2？为何不将（1）、(2)两式写成如下的(3)、(4)两式：
                （3）
                                                           （4）
如果按(3)、(4)两式计算，则模拟的谐波项最低、最高频率将分别为f_min，f_max，这样不是更方便一些而且可以使模拟的随机过程的自功率谱密度更加符合给定的自功率谱密度吗？在星谷胜的书中有证明三角级数法为各态历经高斯过程且自功率谱密度函数在N充分大时逼近给定的的自功率谱密度函数的详细证明过程，但我在其中发现f_n和Δf的定义似乎并不影响证明三角级数法的上述性质。我根据上面的(3)、(4)两式也编了程序试验生成限带白噪声，发现用这两式时依然能够生成满足前述性质的随机过程。

请问各位，为何三角级数法要用(1)、(2)式而不用(3)、(4)式呢？

附matlab代码：
function [Sig,Det_f,fs] = GenerateRandomSigSY(NSig,N,t,fl,fn,PSDFunc)
%本函数用三角级数法产生各态历经的高斯平稳随机过程
%NSig是要生成的信号的组数
%N是三角级数的项数
%t是时间列向量
%fl,fn分别是有显著值的频率范围,fl为下限，fn为上限
%PSDFunc是所要模拟的过程的谱密度函数句柄

%Sig是生成的信号
%Det_f是生成的信号的频率间隔（FFT时应有的频率分辨率）
%fs是生成信号的采样频率

    if nargin < nargin(@GenerateRandomSigSY)
        PSDFunc = @(x) (2 * ones(size(x)));
        if nargin < nargin(@GenerateRandomSigSY) - 1
            fn = 25;
            if nargin < nargin(@GenerateRandomSigSY) - 2
                fl = 5;
                if nargin < nargin(@GenerateRandomSigSY) - 3
                    fs = 102.4;
                    T = 10;
                    t = linspace(0,T,fs * T + 1)';
                    t(end) = [];
                    if nargin < nargin(@GenerateRandomSigSY) - 4
                        N = 200;
                        clc; close all;
                        if nargin < nargin(@GenerateRandomSigSY) - 5
                            NSig = 800;
                        end
                    end
                end
            end
        end
    end
   
    fs = 1 / (t(2) - t(1));
    LSig = length(t);                                                        %一组信号的长度
    USig = unifrnd(0,2 * pi,N,NSig);                                        %生成均匀分布的随机数
%    Det_f = (fn - fl) / ;                                                    %按星谷胜书中的方法计算的频率间隔
    Det_f = (fn - fl) / (N - 1);                                            %频率间隔
    k = [1:1:N]';                                                            %项序列
%    fk = fl + (k - 0.5) * Det_f;                                            %按星谷胜书中的方法计算的频率序列
    fk = fl + (k - 1) * Det_f;                                                %频率序列
    Ak = 2 * sqrt(PSDFunc(fk) * Det_f);                                        %幅值序列
    %现在来构造一些矩阵，以便加速后续计算
    Mx_wktn = 2 * pi * fk * t';                                                %构造2*pi*fk*tn矩阵
    Mx_Ak = repmat(Ak,1,LSig);                                                %构造Ak矩阵
   
    Sig = zeros(LSig,NSig);
    parfor n = 1:NSig
        Mx_U = repmat(USig(:,n),1,LSig);                                    %构造均匀分布矩阵
        Sig(:,n) = sum(Mx_Ak .* cos(Mx_wktn + Mx_U),1)';
    end

    save('Sigs.mat','Sig','t','Det_f');
end

kkkttt

个人理解从数学上讲（3）（4）取法并不是不可以
但是从共振原理的角度来讲，取主控频段的中间值更为合理，因此一般情况下都采用（1）（2）的方法

hcharlie

(1),(2)法和(3),(4)法没有本质的区别，它们生成的都叫伪随机系列，不是真随机，真随机还要在此基础上加时域随机化处理，从有限数列生成无限的真随机序列。
从数列生成方法的简易来看，(3),(4)法更容易实现，从功率谱密度变成频谱，加随机相位，做IFFT即得，速度快多了。
举例：
100~200Hz，N=100，(3),(4)法，从100Hz算到199Hz；
相当于用(1),(2)法时的99.5Hz~199.5Hz的结果，也是从100算到199。
能量移了半条谱线0.5Hz，(3)(4)法简单得多。