路过强答一番。
就举一个例子。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pink Noise Generation with MATLAB Implementation %
% %
% Author: M.Sc. Eng. Hristo Zhivomirov 07/30/13 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = pinknoise(N)
% function: y = pinknoise(N)
% N - number of samples to be returned in row vector
% y - row vector of pink (flicker) noise samples
% The function generates a sequence of pink (flicker) noise samples.
% Pink noise has equal energy in all octaves (or similar log bundles) of frequency.
% In terms of power at a constant bandwidth, pink noise falls off at 3 dB per octave.
% define the length of the vector
% ensure that the M is even
if rem(N,2)
M = N+1;
else
M = N;
end
% generate white noise
x = randn(1, M);
% FFT
X = fft(x);
% prepare a vector for 1/f multiplication
NumUniquePts = M/2 + 1;
n = 1:NumUniquePts;
n = sqrt(n);
% multiply the left half of the spectrum so the power spectral density
% is proportional to the frequency by factor 1/f, i.e. the
% amplitudes are proportional to 1/sqrt(f)
X(1:NumUniquePts) = X(1:NumUniquePts)./n;
% prepare a right half of the spectrum - a copy of the left one,
% except the DC component and Nyquist frequency - they are unique
X(NumUniquePts+1:M) = real(X(M/2:-1:2)) -1i*imag(X(M/2:-1:2));
% IFFT
y = ifft(X);
% prepare output vector y
y = real(y(1, 1:N));
% ensure unity standard deviation and zero mean value
y = y - mean(y);
yrms = sqrt(mean(y.^2));
y = y/yrms;
end
https://raw.githubusercontent.com/cortex-lab/MATLAB-tools/master/pinknoise.m
Prof. Hristo Zhivomirov
简介
Professor Zhivomirov is a member of IEEE, Federation of Scientific-Technical Unions in Bulgaria and Union of Scientists in Bulgaria. His research interests include the field of signal processing, electrical and electronics measurements and Matlab. Department of Theory of Electrical Engineering and Measurements in Technical University of Varna
MIT在线教程
https://www.youtube.com/watch?v=iTMn0Kt18tg
【待续】
。。
我来写一个大尺度的吧。
别误会,是真正的大尺度,宇宙尺度的傅里叶变换。
根据宇宙学原理,我们的宇宙是各向同性的,也就是说,从各个方向看过去应该是一样的。但是如果宇宙是完全均匀各向同性的话,就不会产生我们现在观测到的各种大尺度结构了。我们都知道是随机涨落产生了现在宇宙中各种自引力束缚结构的种子。但是,除了随机涨落,宇宙的密度涨落还有别的什么成分吗?因此,天文学家们就想着通过测量宇宙的密度涨落来研究宇宙的结构和演化,并且成功发现了著名的重子声波振荡(BAO):
那么,要怎么计算密度涨落呢?答案是:数星系。
对,你没看错,就是这么简单。星系作为宇宙中可见物质的重要零部件,很好的示踪了宇宙中的密度分布。因此,星系中的两点相关函数,也就被广泛用于计算宇宙密度场。所谓的两点相关函数,就是在距离为r处发现两个星系的概率,公式长这样:
就是密度涨落场
那么,傅里叶变换在里面起到什么作用呢?这就要提到功率谱了。傅里叶变换是把时域的信号转换成频域的信号,功率谱就是频域信号振幅的平方:
其中的 是 的傅里叶变换:
是某个比较大的体积,在其中 是周期性的。
由于采样的星系是离散的,要转换成连续分布的密度场还需要引入一个泊松误差1/N(shot noise):
的上标是离散(discrete)的意思。N=n ,n是平均数密度。其中 长这样:
是 体积内的星系数目,由于 取得非常小, 只有0,1两个取值。
用这个式子就能算宇宙密度场的功率谱了。
当然,由于实际操作中星系数量十分巨大,直接用上面的式子计算是不太现实的,所以需要用到快速傅里叶变换(FFT),把巡天区域分成离散的网格点,然后用某种窗口函数W做weighting,把离散的星系点做平滑之后再用FFT计算功率谱。
这是另一个故事了,我也还在学习中,什么时候搞懂了再写吧。
最后,这个视频是我见过的最通俗易懂的傅里叶变换的讲解了,这个up的其他作品都很好,墙裂推荐!!!
这个问题简直是为了恒星震动而生的!
恒星是一团很热的等离子气体,在引力和气体压力的平衡下保持稳定。但是,如果其内部有一些扰动,那么恒星的结构就会发生周期性的震动,这就是恒星震动。
恒星震动的一个结果就是,恒星的亮度会发生周期性的变化。所以天文学家通过观测恒星亮度的变化,即可知道恒星是否在震动,也可以把观测和理论做对比,探测恒星内部结构。
现在以一颗星为例,演示一下傅里叶变换在星震里的作用。
开普勒卫星很厉害,对着天鹅座和天琴座之间的一片天区连续观测了四年,天文学家下载好了某一颗星的亮度数据,也就是光变曲线,如下图
横轴为时间,纵轴为亮度。但啥也看不出来是不是?那我们放大看看
放大看好像能看出一点问题了,似乎恒星的亮度在上下波动?但这个波动是真的吗?还是只是测量的不确定度造成的呢?于是我们做个傅里叶变换。
具体怎么做傅里叶变换可以看这篇文章。
做完傅里叶变换后,我们终于看到了一些东西。。
上图就是把光变曲线做完傅里叶变换后的样子,即功率谱。其横轴是频率,单位是天分之一,纵轴是振幅。简单来说,如果光变曲线里有一个频率为 的信号,那么在功率谱里频率 的地方就会有一个峰。
功率谱明显分成两部分,左边的低频信号和右边的高频信号。低频信号(在1天分之一)是恒星内核的震动,回复力为浮力,我们叫它们重力模式。高频信号(18天分之一那里)是恒星外层的震动,回复力是压力,我们叫它们压力模式。
我们先来看压力模式:
上图展示了压力模式某一个震动频率的分裂。分裂的原因是表面的自转。由于自转的原因,震动频率会分裂成三个,分裂大小即等于表面自转频率。所以这颗星的表面的自转大约就是0.01天分之一,即大约100天转一圈。
现在再说说引力模式。引力模式发生在恒星内部,其震动频率较低。上图即为这颗星的引力模式,可以看到每一个峰都分裂成了两个,而且有一系列的峰。分裂成两个即为内核自转造成的,而所谓一系列的峰为不同径向节点数的泛音。。好吧我也不知道怎么举例了。。
理论上预言,引力模式一系列峰,以周期为单位的话,应该会等间隔。我们画一下看看
上图还是功率谱,但是轴坐标为周期了,即频率的倒数。现在可以看到,引力模式的一系列峰确实是近似等间隔的,都在2100秒上下波动,有波动的原因是内部元素梯度。所以天文学家就可以用这个来推算内部氢燃烧阶段,以及元素混合的程度。
自转分裂用加号和减号标注出来了。对于引力模式,自转分裂等于自转频率的一半,算出来这颗恒星的内核自转速度也差不多是100天一圈,和表面一样。
所以我们用压力模式简单算了下表面的自转,又用引力模式算了内核的自转,发现都是一百天。这颗恒星这么巨大,也都是气体不是固体,自转速度从里向外居然是一样的!
这个还是个功率谱,还是引力模式,但是你会发现这一系列的峰(红色点点)的间隔变得越来越小了,这是因为内核自转太快了,大约是0.8天。
上图是我们太阳的功率谱,你会发现这些峰在频率上几乎是等间隔,而且会一组一组地出现。
恒星不是不变的,恒星也会震动。但是直接看震动信号的话,很难看出什么,做过傅里叶变换后,就可以看到一个个的震动频率了。天文学家就可以用这些震动频率去推测恒星内部的环境,也可以检测当前物理模型是不是正确。