Daily Archives: 2013年4月23日

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.2
randn(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_cut
f/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((2
i-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