Gxx1, w = periodogram(vq; nargout = 2);

Gxx1, w = periodogram(vq, rectwin(N), N, fs, "onesided", nargout=2);

using TyMath
using TyBase
using TyPlot
using TySignalProcessing
clear()
clc()
Gqn0 = 1024;
#D级功率谱10^(-6)m^3A=16B=64C=256D=1024E=4096
n0 = 0.1;
#参考空间频率1/m
u = 12.5;
#车速m/s12.5*3.6=45Km/h
Gdq = 4 * pi^2 * Gqn0 * 10^(-6) * n0^2 * u;#速度功率谱
################################################################
fs = 1000;
#单位：HZ
T = 100;
#单位：S
fi = 0.1;
#最小截止频率
fa = 20;
#最大截止频率
L = T * fs;
#信号点数
df = 1 / T;
freq = 0:df:(fs/2);
ni = round(Int, fi / df + 1);#四舍五入取整求最小截止频率对应数组元素的下标
na = round(Int, fa / df + 1);#四舍五入取整求最大截止频率对应数组元素的下标
Gx = zeros(1, Int(L / 2) + 1);#建立一个长度为nft/2元素全为0的向量
Gx[ni:na] = Gdq * ones(1, na - ni + 1);#将r的正频率带通内的元素赋值，生成赋值谱
figure()
subplot(311);
loglog(freq, Gx');
grid("on");
xlabel("频率（Hz)");
ylabel("速度功率谱（m^2/s)");
#################################################################
L = length(freq);
#正频率长度，是信号长度一半+1
N = (L - 1) * 2;
#信号长度
Gxx = Gx;
#双边功率谱密度
Gxx[2:L-1] = Gx[2:L-1] / 2;
#把单边功率谱密度幅值变为双边功率谱密度幅值
Vx = sqrt.(Gxx * N * fs);
#计算期望双边频谱幅值
#用随机数构成相角
fik = pi * randn(1, L);
#产生随机相位角
fik[1] = 0;
fik[end] = 0;
#产生单边频谱(使第一个与最后一个为实频谱逆变换后才能变为实数）
Vk = Vx .* exp.(1im * fik);
#产生单边频谱conj.(Vk[2:L-1)
Vk1 = reverse(Vk);
Vk2 = Vk1[2:L-1]';
#提取后转化为列向量，因此需要（共轭）转
VVk = [Vk Vk2];
#利用共轭对称性求出双边频谱
vq = ifft(VVk);
#傅里叶逆变换
t = (0:(N-1)) / fs;
#时间序列
vq = real.(vq);
subplot(312);
plot(t, vq);
grid("on");
xlabel("时间(s)");
ylabel("路面激励速度（m/s)");
filename = "Exc.jld2";
save(filename; vq, t);
dw = 2 * pi * df;#计算圆频率间隔(rad/s)
w0 = dw:dw:(2*pi*fs/2+dw+dw);
Dk = Vk ./ (1im * w0');
Dk1 = reverse(Dk);
Dk2 = Dk1[2:L-1]';
#提取后转化为列向量，因此需要（共轭）转置成行向量
DDk = [Dk Dk2];
#利用共轭对称性求出双边频谱
xq = ifft(DDk);
#傅里叶逆变换

xq = real.(xq);
subplot(313);
plot(t, xq);
grid("on");
xlabel("时间（s)");
ylabel("路面激励位移（m)");
filename = "Exc_xvt.jld2";
save(filename; xq, vq, t);
#对随机序列x求周期图法的功率谱密度
Gxx1, w = periodogram(vq, rectwin(N), N, fs, "onesided", nargout=2);
fw = w * fs / (2 * pi);
#从[pi,pi]到-fs/2,fs/2]
Gx1 = 2 * Gxx1 * 2 * pi * 0.001;
#先转化为单边谱，再转化为频率
figure()
plot(fw, Gx1);
grid("on");
xlabel("频率（Hz)");
ylabel("速度功率谱（m^2/s)");