Monthly Archives: April 2013

Fourier变换（二）

[mathjax]
上一篇介绍了用Fourier变换来重构一个未知正弦信号的方法。本篇介绍对于未知非正弦信号的情况。

大幅振荡剪切的参数

大幅振荡剪切实验产生的波形是非正弦信号，非正弦是非线性粘弹性的特点。所以对非线性粘弹性程度的表征，就是对波形正弦形状偏离程度的表征。关于大幅振荡剪切和非线性粘弹性的原理与实践，可见Hyhn等的综述[1]，或者参阅我的博士学位论文第1.2.1节、第1.3节和第2章。本文只解决MATLAB编程实现问题。

大幅振荡剪切（LAOS）实验中，我们从流变仪同时获取两个信号。一个是应变的信号，一个是应力的信号（事实上，一个是马达位移信号，另一个是力矩信号，它们都要乘以一个系数，才是应变和应力。此文先略去这一事实了）。在应变控制的实验中，我们施加一个正弦的应变\gamma\left(t\right)，测量应力的响应\sigma\left(t\right)，即
$$\begin{cases}\gamma\left(t\right)=\sin\omega t\\ \sigma\left(t\right)=\sum_{n=0:\text{odd}}^\infty\sigma_n\sin\left(\omega t+\delta_n\right)\end{cases}$$由此式可定义LAOS实验结果的表征参数I_{n/1}\equiv\frac{\sigma_n}{\sigma_1}，以及\mathit{\Phi}_n=\delta_n-n\delta_1，来表征非线性粘弹性的程度（即\sigma\left(t\right)偏离正弦函数的程度。此外，我们还会以\sigma\left(t\right)对\gamma\left(t\right)作图（参数方程作图），得到Lissajous曲线。然后用曲线的形状特征来表征非线性粘弹性的程度（即曲线偏离椭圆的程度）。描述曲线形状特征，往往就是各种各样的切线斜率（可自行上网搜索复习参数方程求导）。在LAOS中常用：$$G_\text{M}=\left.\frac{d\sigma}{d\gamma}\right|_{\gamma=0}$$$$G_\text{L}=\frac{\sigma{\gamma_0}}{\gamma_0}$$$$G_\text{k}=\left.\frac{d\sigma}{d\gamma}\right|_{\sigma=0}$$

MATLAB编程步骤

在详细介绍MATLAB代码之前，先总结一下我们要完成什么任务，跟从哪些步骤。我们的任务就是要从应力和应变的信号计算出上述的这么多参数。实验时，我们会规定应变的幅度\gamma_0和频率\omega，但是，由于仪器误差，实际施加在样品上的应变幅度和频率，未必完全等于的命令值。尽管这一误差非常小，但从做法上我们还是要以实际应变信号为准。所以，我们首先要对实际应变信号进行判断，获得实际的应变幅度和频率。如何判断一个未知正弦信号的振幅、频率和相位的方法已经在上一篇文章中介绍了。

然后，就要处理应力的信号。为了获得\sigma_n和\delta_n同样需要对\sigma\left(t\right)进行fft处理，方法跟上一篇文章介绍的类似，区别只在于我们除了要获得基频的幅值和相位之外，还要获得高次谐波的幅值和相位。理论上来说n=\infty，但是实际上我们根据需要取有限的n个高次谐波，就可以比较准确地重构出\sigma\left(t\right)的波形了。提醒一下，定义式里应变是一个零相位的标准正弦函数，但实际重构的应变波形可能会有非零的相位，这时，应力的各个谐波的相位要都减去应变的相位，才是我们要求的\delta_n。

Lissajous曲线的各种导数，要同时在\gamma\left(t\right)和\sigma\left(t\right)都重构了之后，再进行求导计算得到。总而言之，要完成以下三个步骤：

  1. 重构应变\gamma\left(t\right)，获得\gamma_0、\omega和相位
  2. 重构\sigma\left(t\right)，获得\sigma_n和\delta_n，计算I_{n/1}
  3. 用重构的\sigma\left(t\right)和\gamma\left(t\right)作Lissajous曲线，以及计算G_\text{M}、G_\text{L}和G_\text{k}

重构\gamma\left(t\right)和\sigma\left(t\right)

MATLAB帮助文件里fft命令的条目有一个范例，应该能解决大部分问题了。其中在coherent sampling（见上篇文章）处，范例用next power of 2，其实没有必要。下面的代码，我就不多作解释了。


clear all
clc

% the strain signal
fs=1000; % sampling frequency
t=0:1/fs:30.2; % the time series of the strain signal
t=t(:); % convert to colume vector

w=5; % (rad/s) the angular frequency of the strain signal
f=w/(2*pi); % (Hz)

x=250sin(2pift+1)+25randn(size(t)); % strain signal with noise
y=1.001+8.1532.sin(2pif.t+2.5936)+2.4952.sin(6pif.t+0.5199)+...
0.8856.sin(10pif.t+4.3353)+0.2932.sin(14pif.t+1.8641)+...
0.0992.sin(18pif.t-0.6417)+0.0282sin(22pif.+3.3818)+...
0.2randn(size(t)); % stress signal with noise

[AX,H1,H2]=plotyy(t,x,t,y,'plot'); % Figure 1
set(AX(1),'XLim',[0 8])
set(AX(2),'XLim',[0 8]) % set the time axes

Figure 1 Parts of the strain and the stress signals

Figure 1 Parts of the strain and the stress signals

两个信号都加入了噪音，而且应力信号加入了直流分量，更接近真实的情况。另一个真实的情况是，应力响应的波形往往需要经过几个周期之后才达到稳定的形状。我们当然是取稳定之后的波形来进行Fourier变换，因此就要截取后几个周期。


N=length(x);
N_per_cycle = fs/f; % average points per cycle
N_cycles = floor(N / N_per_cycle); % maximum number of complete cycles


N_cut=round(15*N_per_cycle); % use the last 15 cycles
t=t(length(t)-N_cut:end);
x=x(length(x)-N_cut:end);
y=y(length(y)-N_cut:end);
% crop the last N_cut data.

接下来，先对应变信号进行重构，方法跟上篇文章一样。


Xraw=2fft(x)/N_cut;
Xraw=Xraw(1:ceil(N_cut/2));
freq=0:fs/N_cut:fs/2-fs/N_cut;
A = abs(Xraw(round(N_cutf/fs)+1))
delta = wrapTo2Pi(angle(Xraw(round(N_cutf/fs)+1))+0.5pi)
frq = freq(round(N_cutf/fs)+1)
x_rec=A.sin(2pifrq.*t+delta);
plot(t,x,'b-',t,x_rec,'-r')

Figure 2 Reconstruct the strain signal

Figure 2 Reconstruct the strain signal

然后轮到应力的波形。先看看它的Fourier变换结果：


Yraw=2*fft(y)/N_cut;
Yraw=Yraw(1:ceil(N_cut/2));


plot(freq, abs(Yraw)) % Figure 3

Figure 3 Amplitude spectrum of the stress signal

Figure 3 Amplitude spectrum of the stress signal

放大最左边的部分，可以看到好多高次谐波。为什么右边有这么长的多余部分？回忆上篇文章。Fourier变换后的频率范围等于fs/2。我们fs达到1000 Hz，而信号主频只有约0.8 Hz，几个高次谐波都集中在最左侧了。有时间的话，最好多点尝试不同信号的FFT，熟悉这一转换是认识大自然的一个新的思想高度。因为大自然很多自然过程就是在做Fourier变换。例如X线穿过晶体发生衍射，光就被Fourier变换了，形成的二维斑点图可通过进行二维的反Fourier变换，还原成样品中的晶格排布。所以我们总是说，X射线的图谱是“倒易空间”（Fourire变换前后坐标变倒数了）。

Figure 3还可看到，在频率为0处有一个峰。这个峰是原波形的直流分量（见文章开头构建信号的式子），我们需要知道各个谐波在信号序列中的具体位置。上篇文章已经解决了正弦信号的情况，对于高次谐波（都是奇数倍频），只要相应地乘以3、5、7即可。但是，根据定义式，应力的Fourier展开理论上是无穷多项的，而实际上也需要取尽量多项以确保最准确地重构原波形。到底Fourier变换之后，最多允许看到多少个高次谐波呢？这就相当于要找使nf<fs/2的最大奇数n。不过，如果信号的取样频率fs相当大，最大谐波次数n也会非常大，我们实在没必要取这么大的n（如Figure 3的情况右边大部分是多余的。因此我们又人为规定最大谐波次数不超过25个，就是说我们假定25个高次项对于再怎么非线性的信号都肯定足够了。以下是相关的代码：


max_nr_harm=floor(fs/(2*f)); % Maximum number of harmonics
if max_nr_harm/2==floor(max_nr_harm/2)
max_nr_harm=max_nr_harm-1;
end
% maxmum odd number of harmonic

max_nr_harm=max_nr_harm-2; % fall back one more odd harmonic for safety

if max_nr_harm>25
max_nr_harm=25;
end

接下来就是要确定这些谐波的幅值和辐角：


% create the variables for harmonic amplitudes and phase angles
I_n=zeros((max_nr_harm+1)/2,1);
delta_n=zeros((max_nr_harm+1)/2,1);


for i=1:(max_nr_harm+1)/2
I_n(i)=abs(Yraw(round((2i-1)N_cutf/fs)+1));
delta_n(i)=angle(Yraw(round((2i-1)N_cutf/fs)+1))+0.5pi;
end
DC=abs(Yraw(1))cos(angle(Yraw(1)))/2;

其中，直流分量的大小计算方法为什么是这个，我不太记得了。这时，你已经可以通过I_n(i)和delta_n(i)来计算相应的非线性粘弹性参数I_{n/1}\equiv\frac{\sigma_n}{\sigma_1}和\mathit{\Phi}_n=\delta_n-n\delta_1了。再次提醒，应力波形的各谐流相位角要减去应变波形的相位角，才是定义式中的相位角。我们也可以检验一下重构的应力波形：


y_rec=zeros(length(y),1); % create the variables for the reconstructed signals
for i=1:2:max_nr_harm
y_rec=y_rec+I_n((i+1)/2).sin(i2pif.*t+delta_n((i+1)/2));
end
y_rec=y_rec+DC;
plot(t,y,'-b',t,y_rec,'-r') % Figure 4

Figure 4 Reconstructed stress signal

Figure 4 Reconstructed stress signal

作成Lissajous曲线来看：


plot(x,y,'b-',x_rec,y_rec,'r-')

Figure 5 Reconstructed Lissajous curve

Figure 5 Reconstructed Lissajous curve

可见，通过Fourier变换重构波形，可以很好地把噪音排除掉。Lissajous曲线的各个形状参数，需要在曲线上找各种交点求切线（见定义式）。如果在原信号的基础上直接做，可想而知误差会有多大。这就是为什么我们明明已经有信号不直接用要先重构一番的原因。具体如何在Lissajous曲线上找点、求切线的方法，在下一篇文章中介绍。

References

  1. K. Hyun, M. Wilhelm, C.O. Klein, K.S. Cho, J.G. Nam, K.H. Ahn, S.J. Lee, R.H. Ewoldt, and G.H. McKinley, "A review of nonlinear oscillatory shear tests: Analysis and application of large amplitude oscillatory shear (LAOS)", Progress in Polymer Science, vol. 36, pp. 1697-1753, 2011. http://dx.doi.org/10.1016/j.progpolymsci.2011.02.002

知乎：对于搞科研学术的研究生而言，究竟该专注于平心静气做好研究工作，还是应该为以后找工作早做打算？

原问题地址：http://www.zhihu.com/question/20949946

看了本面其他人的回答，我想补充两点。反正我的回答已经长得不行了，再长点也就那样。

一、我相信这个问题答案不应该是直接让提问者做哪个选择。再加上，我在知乎上回答问题都是借问题来表达一个观点——即那种既适用于具体这个问题，又有点普适性的道理。因为，我希望自己是一个讲原则的人（或者说读书读傻了人死板），我做这个选择的理由，跟我做那个选择的理由，要能够同时符合一个自洽的价值观体系，而不是见利忘义，此一时彼一时，见人说人话见鬼说鬼话……如此各般等等。不是说现在这么做对我有利所有我这么做，而是我一向信奉什么所以我这么做。成功学畅销书也许能帮你分析利弊，但不能强加你一个哲学。人生不是一场生意，不要去做人生的商人。

做同样一个选择，可以基于很多理由。我不赞同你做这个选择的理由，不等于我不赞同你做这个选择。我批评那具理由，不等于我批评那个选择。若要问“那你到底认为该做何种选择呢？”我当然不说，也不应该说，选择就是应该自己做的。我的任务就是摆一通理由，而不是直接说，你就选A吧B是错的。

读了这么多书，至少要形成“知其然知其所以然”的习惯，而不是还停留在“肚子饿了要吃饭，困了要睡觉，来性子了就撸一撸，地上有钱了就捡，枪指头上喊救命”的本能反应水平。只有未成年人才公认不具备承担责任的能力。成年人就是要为自己的选择负责。

在硕士和博士阶段认识“科研是什么”，完全不晚。我想不到有什么方式可以更早地认识这个问题。只是人完全可以选择压根不去认识这个问题。既然你对这个问题产生好奇，你去认识了一番，那就不应该又去抱怨问题的答案让你失望，更不应该说是什么“社会上专家学者”忽悠了你。我天天说哪个歌手很好听，你觉得我从小弹钢琴懂音乐，于是你自己买了张演唱会票，回来你能怪我么？只要我没直接叫你去买演唱会票，我又没有收那歌手的钱违心地宣传他（真心觉得好听），我就没有过错。这不能说是我狡猾，也不能说你被忽悠了，而完全是个正常的代价。

二、我的答案的重点在于“要和先导师商量、讨论”。这也是匿名用户建议的第一条。主要是由于我估计跟导师讨论之后的结果多半就是只能维持做实验的强度，我的答案才会显得好像在叫他将来选择做科研。但我不排除你和导师商量了之后导师就放你跑了这种可能。匿名用户的结局就是这么好。至少，你先跟导师说一声，导师会有所准备，假如在课题项目上安排得过来，说不定也可能允许你缺勤。匿名用户所谓的：“和你的导师谈。他能接受你不做学术，你的朋友有什么好质疑你的？”换个说法就是，只要你有本事说服你导师接受你不做学术就万事OK；即如果说服不了的话就只能自己兜着走。

这跟你不先商量自己跑了是两码事。我十分想强调的是诚信。可是我发现很多中国人的心态就是，如果明知跟你商量了结果对我不利，我就还不如不跟你商量（哪怕道义上我不应该）。即诚信如果对我有利我就讲诚信，如果对我不利我就不讲了。我首先不赞同这种做法。正确的做法是：明知跟你讲了八成对我不利，不跟你讲我就自己完成了；但我还要跟你讲，然后就真的吃亏。这就叫做“不要见利忘义”。

现在我谈谈这里面有什么“义”。我刚才已经说了，并没有谁忽悠你来读研。既然是研究生，那么就不仅是学校的学生，还是课题组这一团队的成员，按照学校和课题组的规定来做事。例如每天八点半到实验室之类。跟导师商量好的实验计划，导师当然有理由定期追问，你当然被假设是所有的工作时间都在课题事务上，因为你读的是全日制研究生。当初要进课题组让你进了，结果你实验不做，实验室不去，要找你你去了别的城市，明天才能回来，突然连续一个月失踪，原来参加了企业实习生……我认为这不能拿你的前途问题来当借口。你嚷嚷着说科研就是个大忽悠，你被骗了，你根本不喜欢科研，你后悔了；但这些都是在你读研之后，体验过了解过之后才说的话。假如这些是你读研之前就知道的事，那你就更恶劣，你故意假装想读研而考入本课题组，主观上就故意来混学历么？你可以选择将来工作，但就读期间你只能完成导师和你讨论商定的实验计划和其他科研任务（写论文投稿）。而且我特别要说的是在这方面你跟导师讨论的余地不多，你没有太多的理由要求导师因为你将来不做科研而让你现在就消极怠工。

在很多现在的研究生看来，这简直要立马上论坛大骂这个导师了。在他们眼中，读研时做不做实验天然就是可做可不做。因为不同的导师不一样严，显得好像也没有一个准则到底要多努力才算合格。自己要多干事，无非是这一个导师特别严。自己到底是多花点时间在实验上呢，还是不做实验去搞找工作的事呢，完全看自己的选择（于是有了本问题），而不用考虑别人。其实不管严也好，松也好，研究生是课题组的一员，是为完成课题组的研究目标而来的。哪怕所谓“自由探索”，也只能是渗透在课题实施的各种环节，不是说整个课题就是“自由探索”。因此，研究生是完成导师所立科研项目的主力。这并不是“免费”或“廉价”劳动力，而是一个完全不惧公开的“行规”。科研经费申请书形式上是要求写明参与的研究生姓名、身份证号码，一年工作时间（多少个月）等信息的。也就是说公认了课题项目是导师立项研究生去完成的，研究生有这个义务。相应的科研酬金以学生三助的形式发放给研究生。但由于项目和经费不同导师不一样，我们学校的做法是学校统一从财政出一部分，然后导师从经费中出一部分，保证每个研究生的补助金额一样。没经费的导师于是就不能招生。这只是避免不同学生补助金额依赖实际参与课题的起止时间的技术性做法。但道理上是因为你参与了课题你才拿到补助。研究生在课题组还有别的义务，例如维护常用仪器，实验室安全，药品管理等等。你突然有一天后悔了，“醒悟”了，然后就不管不顾了，这对整个课题组都是不负责任的。

骂中国的高等教育是一回事，骂社会上的“专家学者”，骂奇葩导师、骂四六级、骂那些你鄙视的书呆子们、那些“万般皆下品唯有读书高”的傻子们，是一回事。“诚信”和“责任心”是另一回事。这些不要说研究生了，从本科生开始就已经很缺乏。本科申报学生科研计划进实验室，来找老师。问他：你每周都哪些时间有空？答曰：周五晚，周六晚和周日下午。课这么满吗？哦，我XXX时间做家教、XXX是预备党员组织生活、XXX协会每周X例会，XX又是家教。我说，你每周忙到只剩这点时间了，还不休息？来实验室，时间也不够啊。他说：参加这个加很多份，而且不用做实验，只要跟着师兄的论文挂个名，将来能评奖……大四的毕业设计，凡是不读研的，基本为零。老师都只能理解一下因为这也是现在的大环境逼的，但这也是做坏了规矩，恶性循环。有的老师要求严一点，学生拍桌子翻脸。其心态就是，凡是我有什么需求要实现，世界上就不应该有别的事可以挡着我。现在我要考公务员，这才是最王道的。其实，这样的学生，我还真希望快点去工作几年，千万别继续留在学校读研。早几年认识到什么才是王道更好，要是一直留在学校，导师又很好说话，那说不定还真的是读到博士毕业30岁了还是这副德性，人生真毁了。

以下是原来的回答。我的回答中有一句话：要想就业好，就要就业早。不知道为什么，似乎还是给人感觉我是在“为学术多说了几句”，甚至是建议提问者将来去做科研。翻了一下我其他问题的答案，发现，叫年轻人早点接触职场早点实习的观点我不止一次表达，例如如何让注意力集中在单独一件事上？但是，知其然要知其所以然，重点是理由。是理由让我跟别人的观点不同，我之所以为我。
====================================
我觉得匿名用户的答案比我好，不想看长篇大论的，看了他的答案基本就了解我的观点了。

一个人怎么取舍，除了利弊的衡量之外，还有个人性格气质方面的因素。哪怕明摆着有利的选择，你就是不喜欢，也可能不做。你因为自己不喜欢而不做这个有利选择，很可能事后也不会后悔。思考的方法大概就是先客观地分析利弊，再主观地掂量一下。后者是很容易的，自己还不知道自己么，但前者对于阅历尚浅的学生来说往往很难——他太不知道人家。

为以后工作“早”作打算的“早”，是有多早？你读书读了十几年，最后两年才为工作做打算，算不算早？为以后工作早作“打算”的“打算”，又是什么“打算”？好吧，现在你可以不做实验了。那你会做什么呢？每天做的什么事情，才属于“为以后工作早作打算”呢？现在你马上要做什么事情与做实验冲突了吗？毕业工作不好找，你是已经想到什么事情，现在去做，你毕业工作就会好找吗？

现在学生普遍存在的问题，一个是高谈阔论（这其实不仅限于“现在的”学生，而是年轻人普遍喜欢）。二是只做能立马兑现成果的选择，说“功利”已经不够准确了。问题是，这两个问题之间就是矛盾的。真要追求做每一个选择都立马兑现，最要不得的就是高谈阔论，而是要变得非常实际和通透才行。如果你的水平只限于高谈阔论，你休想做什么选择能立马兑现，乖乖地走弯路碰钉子去吧。

做导师的当然会抱怨，现在的学生全都不是来做科研的，都是为了拿学历，至少缓冲就业。做导师抱怨这个，是因为于是他们就没办法对学生要求太高，课题的质量和档次就提不上来。做导师抱怨这个也是很无奈，说到底也是不能怪学生的。因为他不喜欢的话完全可以不招学生，现实是他也必须招，因为生源的普遍情况就是这样，没学生干活不行。这是买方还是卖方市场的问题。

但学生不是完全没有考虑这个问题的义务的。在读研甚至读大学之前，我们的学生太习惯于老师必然是无私奉献的，以教育学生为自己全部的职责的，所以读研之后大多数会对研究生导师不大满意，跟自己以往对老师的印象反差很大。事实上，读研是一种合作。学生要获得培养，导师要完成课题。于是你问：导师为什么就没有培养学生的义务？本来是有的。但现实是，导师只能把学生培养成为独立的科研工作者，培养成企业主管不是导师能力和责任所在，但大多数学生并不愿意接受这样（科研）的培养，所以只好退到“合作”这一步。如果连“合作”都谈不拢，就干脆别选择这个导师了。

可是有很多人事先没有跟导师谈过这个问题，都按自己的想象行事。最后毕业成问题，才突然撕破脸皮。所以，如果你要为你将来就业“早做打算”，第一件事要做的就是问你导师，你毕业的条件是什么。学校的要求满足了是不是就足够了？导师是否有自己的要求。其实你现在谈都迟了，很可能你会发现这个导师的要求很高，你如果要完成他的要求，还真没有时间做别的事情。说不定你导师会说，你现在这个反应做不出来就别想毕业（你做熟手了，毕业走了，上来的学生又是生手。如果你不把成果做到发表那一步，决不放你走是很有可能的，但也不是绝对）。具体，要你自己跟你导师谈。由于现在谈是迟了的，所以谈出来什么结果，你都只能自己吃不了兜着走，只能说是你现在事后了解一下，不是真正的谈判。总之，你处于学生的状态，毕业是首当其冲的重要问题。没有顺利的毕业，就业什么都是浮云。你早做任何其他打算都属于无用功。

不知道你在什么学校，什么城市，什么专业背景，什么研究方向，你的师兄师姐就业就如何“不理想”。你能排除你自己对“职业”的理解还很幼稚这一因素吗？硕士刚毕业是想要有多理想？我不是想深入到这个话题。而是提醒你，如果要为将来就业早做打算，首先要理解就业本身是怎么一回事。你如何制定你的职业规划？职业规划不是凭空制定的。一开始你需要了解很多别人走过的路，例如刚毕业他什么情况，两年后他怎么样，五年后，十年后他又怎么样；从中你发现一些疑问，就能加深你对你这个专业在社会上对应的行业的理解，从而自然而然地形成你自己的职业规划。这些信息，课题组里应该是能了解到的。你先要认识这个行业，了解你们这个专业出去的人，在这个行业中发挥什么作用，扮演什么角色。只要保持对行业持续加深理解和实践，毕业刚开始差一点影响不大。

我不知道你了解完之后，会有什么感觉。我的经历告诉我，除了英语和学历之外，在学校里基本没有什么可以为将来就业早做的“打算”。想要就业好，就要就业早（但你又想要学历，所以去读个研）。工作之后的事情，有八九成是在学校里很难培养到的。千打算万打算不如早点进入社会。有的专业研究生阶段很闲，会去实习。但化学专业要做实验往往还要通宵，这个是谈不上的。不过，你心理也要平衡：为什么要做那么多实验还要通宵？无非是怕毕不了业。而刚都说了毕业是首当其冲的，这你没办法。既然如此，还不如专心把这学历和学位拿到手。

你是否觉得现在做实验，无非是读研阶段要做的事情，将来一辈子都不会再做化学实验了，加上这么困难，所以你现在也不想好好做了。如果你阶段性的事情你就不好好做，培养这种气质属于“为将来就为早作打算”吗？那么用人单位会把你放在一个永不升职、永不调换、永不超生的岗位去——这种老实人只有长期做一件事才愿意好好做，就让他在这个岗位上做一辈子吧。我敢保证，你毕业后的第一份工作一定是不好的，你肯定要换；我也敢保证，你换的第二份工作也不好，因为那时候你有男朋友（女朋友），他她很远，他她这个他她那个不对付；我也敢保证，你换的第三份工作也不好，因为你小孩上幼儿园在XX区要接受，你妈生病了每周四要去XX医院挂专家门诊你必须开车送她，你工作脱不开身要换……但是，如果你第一份工作做得不好，你甚至连第二份这样能不定期支付给你女友搞点小惊喜的工作都找不到，更不要说你的第三份工作那样虽然地点远但是能养活你妈你小孩你全家的工作了。那个时候很可能你早就忘了你研究生课题那个反应物重结晶的温度了。但你把每个阶段性任务做好的素质正在默默地为你的人生保驾护航。我现在讲的还只是做事的素质的其中一个小小的方面。你还在学生时代都嚷嚷说实验很忙百忙之中抽不出时间来做别的，我看你离好的就业单位也比较远。

一个人是生活在一定的人际关系和体制环境之中的。并不是说因此你就人在江湖身不由己，而是说自由就不是天然获得的，而是要寻求，甚至争取的。坐在哪里又不寻求又不争取，只会觉得所有门都是关着的。现实就是你不能全不管你的课题要求又不去实验室跑去参加EF商务英语课程考公务员注册会计师考驾照做企业实习生，但又不能好像永远都读研一样在一个100%不可能拿诺奖的课题上挖掘你自己都不知道是啥的深度结果其他素质全部停留在应届本科生的水平就没提高过。现实就是你必须宁为瓦全不为玉碎。现实就是被迫要做那种妄图鱼与熊掌兼得努力。想法不能停留在怨天尤人或者高谈阔论的层面。到底你在现在这种阶段，要为将来就业做哪些准备？要占用多少时间？你如何安排你的工作和生活去适应？问题具体化，做决策，执行计划，成果管理，培养一下这些事情，也算是“为将来就业早做打算”吧。

最后讲一句：学生的知识面和眼界太窄，没有资格去做那些能立马兑现的打算——这种欲望是毛头小孩都有的，但能力是军师级别人马才具备的。学生能做的基本上都是那种长远起见的积累，读书啦，参加社团啦，学外语啦什么的。不要去鄙视这些事情，因为你只能做得好这些，乖乖地把这些都做好来。麦子成熟的时候，果实饱满的麦穗是低头弯腰的，直杆的都是空壳。

Fourier变换（一）

[mathjax]

为了带师弟做实验，简要总结一下在MATLAB里进行Fourier变换的一些代码，以供大家参考。

Discrete Fourier transform (DFT): The sequence of N complex numbers x0, …, xN-1 is transformed into an N-periodic sequence of complex numbers according to the DFT formula

X_k=\sum_{n=0}^{N-1}x_n e^{-i2\pi kn/N}\quad k=0,1,\cdots,N-1.

一般把变换前的信号xn当成时域信号，变换后的Xk当成频域信号。可见，时域是N个单元的信号，变换后频域还是N个单元的信号。什么周期啊频率啊这些信息是暂时没有的。例如xn为以下这个序列（Fig. 1），横坐标就是序号n。这段序列，可以是一段10秒钟的的信号，也可以是一段100秒钟的信号。当然，实际情况下一段信号是多少秒是事先知道的，出来的信号一共有多少个单元值（N）则取决于我们隔多长时间取一次样。如果一段10秒钟的信号，我们每秒取一次样，这个序列就有10个单元（N = 10）。反之，Fig. 1的这个序列N = 93，假如这个信号是每隔0.1秒取一次样得到的，这个信号的总时长就是9.2秒；而假如这个信号是每隔10秒取一次样得到的，这个信号总时长就是920秒了。每隔多长时间取一次样，可以用取样频率fs来表示。fs = 10 Hz就是每隔0.1 s取一次样。所以，为了给一对离散Fourier变换xn和Xk赋以周期或频率的信息，必须给时域信号规定一个fs。一旦fs规定了，周期和频率也就规定了。还拿Fig. 1的序列为例，设该序列的取样频率fs = 10 Hz，则该序列每个单元隔了0.1 s。我们可以数一下这个序列一个周期有多少个单元——20个，那么该序列的周期T = 2 s，频率f = 0.5 Hz。

Figure 1

Figure 1

事实上Fig. 1的序列是一个有一定时延的正弦序列（即不是从零开始），我们已经根据规定的fs知道了它的周期和频率，不仿就把它改作在时间轴上（Fig. 2）。

Figure 2

Figure 2

可以看到，该正弦信号的零位在第17和第18个单元之间，即在1.6 s到1.7 s处，以1.6 s计算，该信号可写成x\left(t\right)=A\sin(\left[2\pi f\left(t-1.6\right)\right]=A\sin(\left(2\pi ft+\delta\right)，则相位角\delta=-1.6\pi，换成0到2π之间就是\delta=0.4\pi=1.2566。若取时移是1.7 s得\delta = 0.3\pi = 0.9425，所以δ应该在0.9425到1.2566之间。从图中还能读出振幅A = 2，则信号x\left(t\right)=2\sin\left(\pi t+0.9425\sim1.2566\right)。用正弦函数去拟合Fig. 2的序列发现其实最精确的函数描述应该是x\left(t\right)=2\sin\left(\pi t+1\right)，见Fig. 3。

Figure 3

Figure 3

现在我们试着在MATLAB通过Fourier变换来精确地得到xn的幅值和相位，重构一个xn所对应的连续信号。首先构造上述的时域信号：


t = 0:0.1:9.2; % fs = 10 Hz
t = t(:); % convert into column vector
f = 0.5; % f = 0.5 Hz
x = 2 * sin(2 * pi * fi * t + 1)

前面说过对一组N个单元的序列进行Fourier变换之后得到的Xk也是一个N个复数单元的序列。我们看看它长什么样（分开模序列和辐角序列）：


Xraw = fft(x);
stem(abs(Xraw),'o'); % Figure 4
stem(angle(Xraw),'o'); % Figure 5

Figure 4

Figure 4

Figure 5

Figure 5

我们从Fig. 4可以看到结果跟我们想像的差很远。主要有以下四个问题：

  1. 频域谱分成了对称的两半。要取左边一半。但总单元个数为奇数（N = 93），到底是取46个还是47个？
  2. 若只看一半，频谱出现了单峰是没错（因为我们的信号只有一个频率），但它不是只在一个频率上，而是有一个很宽的分布，似乎我的输入信号的频率是“不纯”的——这并不符合事实
  3. 幅值不对，信号幅值是2，但频谱的峰高快到80了
  4. 横轴如何对应到真实的频率轴？

我们一个一个问题地解决。第一个问题，按照fft算法的规律，若N为偶数，则左右分成各N/2的两半即可；若N为奇数，则处于中央的那个单元是左右两半谱共用的。所以，无论取左半还是右半，都取(N+1)/2+1个。综合奇数和偶数的情况，对任意正整数N的取法就是对N/2向上取整。


N = length(x);
Xraw = Xraw(1:ceil(N/2));
stem(abs(Xraw));

第二个问题，之所以频谱的峰有一个很宽的分布，是由于谱泄漏（spectral leakage）现象，简单地说就是我们用来做Fourier变换的信号序列长度并不是恰好整数个周期。只有截取整数个周期来做Fourier变换，才能避免谱泄漏。可是对于一个未知频率的信号，我们本来就是要通过Fourier变换来知道它的频率，又怎么能在此之前就得知信号的周期呢？这种情况下，只能采用另一种方法，那就是把信号的首尾进行某种衰减处理，而这一处理在数学上又要保证不严重影响Fourier变换的结果，称为加窗操作（windowing）。详细的原理不在此详述。作为示范，我们可以利用我们对时域信号的已知信息来解决。已知信号周期T = 2 s，取样频率fs = 10 Hz，到底截取多少个点算一个周期呢？前面计算过了，20个点。干脆我们就只取一个周期得了，试试看：


x_cut = x(1:20); % N = 20
Xraw = fft(x_cut);
Xraw = Xraw(1:10); % N=20, so ceil(N/2) = 10
stem(abs(Xraw)); % Figure 6

Figure 6

Figure 6

这次，频峰值终于只出现在一个频率上了。而由于用来做Fourier变换的信号N减小到了20个，所以频域信号也就只有20个单元，取左半边就是十个单元。但我们看到，峰值也变了，原来近80，现在是在20。这说明，这个不正确的峰值，是跟N值有关的。这就涉及到第三个问题：峰值怎么校正回来。其实，我们做离散Fourier变换，需要进行归一化。但是怎么归一化，要看离散Fourier变换和反Fourier变换的具体规定，不同的书规定不一样。MATLAB里fft命令是需要用N来归一化的。但如果出来双面的频谱，各边就要用N/2来归一化。所以进行了归一化的的Fourier变换操作应为：


Xraw = 2 * fft(x_cut) / 20;
Xraw = Xraw(1:10);
stem(abs(Xraw));

出来的结果，峰值就会准确落在2（图略）。现在剩下最后一个问题，那就是横轴怎样换成真正的频率？这就要先了解Shannon-Nyquist定理。在此我只把结论说出来：频域信号单元的间格频率是fs/N。因此整个左半边谱的频率范围就是0 Hz ~ fs/2 – fs/N。为了解决奇数偶数问题（即以防fs/2不一定是fs/N的整数倍），最终我们就如下构造相应的频率轴：


freq = 0:fs/20:(ceil(N/2)-1)*fs/N;
freq = freq(:);
stem(freq,abs(Xraw)); % Figure 7

Figure 7

Figure 7

其实，要取整数个周期，不一定要只取一个周期。原时域信号至少有4个周期（见Fig. 1），所以我们可以取满整4个周期来做Fourier变换：


x_cut = x(1:80); % 4 cycles = 80 points
Xraw = 2 * fft(x_cut) / 80; % N = 80
Xraw = Xraw (1:40);
freq = 0:fs/80:fs/2-1/80;
stem(freq,abs(Xraw)); % Figure 8

Figure 8

Figure 8

可见，当我们多取几个周期时，Fourier变换的结果频率“分辨率”更高。因此为了获得尽可能完整地信息，要在保持取整个周期的同时，尽量取最多个周期来进行Fourier变换。我们把相应的辐角谱也作出来：


stem(freq,angle(Xraw); % Figure 9

Figure 9

Figure 9

我们发现，跟Fig. 5相比，Fig. 9的辐角谱显得非常杂乱无章。这是由于我们很好的消除了谱泄漏。规律：谱泄漏越严重，辐角谱越有序；谱泄漏避免得越彻底，辐角谱就越杂乱无章。虽然杂乱无章，但是模谱的幅值所在相应频率对应的辐角是应该是正确的。我们看看0.5 Hz处的辐角：


wrapTo2Pi(angle(Xraw(5))) % wrap result within 0 to 2Pi

结果不对，这是由于在MATLAB的fft得到的辐角差了π/2，要加上：


wrapTo2Pi(angle(Xraw(5))+0.5*pi)

这次就对了。

问题是，在实际应用中，我们不可能做一个Fourier变换，就作一次图，用肉眼看出峰值落在哪个频率上。必须事先就知道信号频率f在频域信号中的对应序号，直接就给出模的峰值和相应的辐角。所以，现在的问题是，如何不通过作图就知道，对于一个频率为f，采样频率为fs，N个单元的信号xn，进行了上述的Fourier变换后得到的Xk，第几个单元呢对应着f？这很简单，Fourier变换后的频率间隔是fs/N，因此，第Nf/fs+1个单元就是模的峰值所在频率。本例子中既可验证。

但是，假如信号的频率f在Fourier变换之后的频率轴上的位置，是恰好在频域信号两个单元之间呢？例如，我的信号频率是f=0.51 Hz，fs和N值不变，按照上述算法，要取第5.08个单元，当然不存在这个非整数单元。这就是说，我们要对Nf/fs+1这个公式加一个四舍五入的取整。注意，由于信号频率不同了，“整数个周期”的截取可能会略有不同。以下是具体的运算过程：


f = 0.51
x = 2 * sin(2 * pi * f * t + 1);
N = length(x);
N_per_cycle = fs/f; % average points per cycle
N_cycles = floor(N / N_per_cycle); % maximum number of complete cycles
N_cut = round(N_cycles * N_per_cycle); % number of points for FFT
x_cut = x(1:N_cut);
Xraw = 2 * fft(x_cut) / N_cut;
Xraw = Xraw(1:ceil(N_cut/2));
freq = 0:fs/N_cut:fs/2-fs/N_cut;
abs(Xraw(round(N_cut*f/fs)+1))
wrapTo2Pi(angle(Xraw(round(N_cut*f/fs)+1))+0.5*pi)
freq(round(N_cut*f/fs)+1)

结果离真实值就稍有偏差。这时，如果我们试试把fs提高到1000 Hz。


clear all;
fs = 1000;
t = 0:1/fs:9.2;
t=t(:);
f = 0.51
x = 2 * sin(2 * pi * f * t + 1);
N = length(x);
N_per_cycle = fs/f; % average points per cycle
N_cycles = floor(N / N_per_cycle); % maximum number of complete cycles
N_cut = round(N_cycles * N_per_cycle); % number of points for FFT
x_cut = x(1:N_cut);
Xraw = 2 * fft(x_cut) / N_cut;
Xraw = Xraw(1:ceil(N_cut/2));
freq = 0:fs/N_cut:fs/2-fs/N_cut;
A = abs(Xraw(round(N_cut*f/fs)+1))
delta = wrapTo2Pi(angle(Xraw(round(N_cut*f/fs)+1))+0.5*pi)
frq = freq(round(N_cut*f/fs)+1)

结果就比之前更接近真实值了。这说明，时域信号的采样频率越高，越有利于精确地重构。现在可以看一下重构的效果：


x_rec = A * sin(2*pi*frq*t+delta);
plot(t,x,'bo',t,x_rec,'r-') % Figure 10

Figure 10

Figure 10

下一篇文章会介绍当原时域信号的频率“不纯”的时候如何通过Fourier变换进行重构。