clear all;
close all;
fd = 10;
%需要2″的時長，則 N/(2*fd)>20
N  = 2^( floor(log2(40*fd))+1 );
%根據頻譜采樣得到的時間間隔為1/fd, 需要1″1000個樣點，則需要時間間隔為1/1000s，因此需要插入點為
f_sample=floor(1000/fd)+1;

%產生高斯隨機變量
x= randn(1,N/2)+j*randn(1,N/2);
y= randn(1,N/2)+j*randn(1,N/2);

%頻率濾波
inx = [conj(x(end:-1:1)) x];
iny = [conj(y(end:-1:1)) y];

df = 2*(fd-0.01)/(N-1);
f = -(fd-0.01):df:(fd-0.01);
fdfilter = 1./sqrt(sqrt(1-(f/fd).^2));%.*exp(-j*2*pi*f);
%fdfilter = fftshift(fdfilter);

lx = inx.*(fdfilter);
ly = iny.*(fdfilter);

zx = [lx(N/2+1:end) zeros(1,f_sample*(N-1)) lx(1:N/2) ];
zy = [ly(N/2+1:end) zeros(1,f_sample*(N-1)) ly(1:N/2) ];;
%zx=lx;
%zy=ly;

cx = real(ifft((zx)));
cy = real(ifft((zy)));

r = sqrt(cx.^2+cy.^2);
%功率歸一化
r2 = mean(r.^2);
r = r./sqrt(r2);
figure(1);
plot(r);
axis([1 2000 0 max(r(1:2000))])