为LAOS波形重构做了一个GUI。
波形重构GUI
除了上一篇文章提到的一些问题外，还有其他的难题。最主要的就是，当原信号的数据量特别大的时候，正弦波拟合函数和MITLaos都会出错。我设计了一个判断，适当时候把信号down sample了再做正弦波拟合和MITLaos，前者基本没问题了，但后者还是不时出问题。只要正弦波拟合没问题的话，波形重构是没问题的，没必要卡在MITLaos上。所以我打算不用MITLaos，直接在重构了的波形上求导得到G_M、G_L等参数。事实上，在程序中G_k参数就是通过直接求导得到的，结果令人满意。
对真实世界的信号进行求导，需要先大幅提高信号的信噪比，否则求导运算就是一个噪音放大器。对于LAOS实验的信号，由于假设其为一系列奇次谐波的和：
。因此，可以据此直接进行波形重构。即通过FFT得到信号的谐波振幅和相位角，按照此式重构一个平滑的波型，再拿去做微分。原理很简单，但是在MATLAB里整了半天，强度是对的，但相位角永远是错的。研究了一晚上，才发现FFT得到的相位角总是少了90°，也不知道什么原因。总之加上这90°之后就完美了。为了减少日后摸索，把代码贴在这里：
clear all;close all;clc;
%============ Generate the signal for reconstruction ============%
fs=57; % sampling frequency
fi0=10/3; % input frequency
t=0:1/fs:123456/fs; % the length of signal is 12345.
t=t';
delta=[1 3 2 1];
x1=sin(2*pi*fi0*t); % the input signal is a sinusoidal wave;
x2=12*sin(2*pi*fi0*t+1) +3*sin(3*2*pi*fi0*t+3)+0.8*sin(5*2*pi*fi0*t+2)+0.02*sin(7*2*pi*fi0*t+1) +2*rand(size(t)); % the output is a sum of odd harmonics
x=[x1 x2];
%============ Reconstruction begins ============%
fi=sinfapm(x(:,1),fs); % the frequency is obtained by sinewave fitting.
Nr_cycles=floor(fi*length(x)/fs); % number of comlete cycles
N=round(Nr_cycles*fs/fi); % cutoff length for coherent sampling
x=x(1:N,:);
t=t(1:N); % cut the time axis too
Xraw=fft(x(:,2));
phase=angle(Xraw);
freq_bin=0:N-1; % frequency bin
freq_res=fs/N; % frequency resolution
freq_i=freq_bin*freq_res; % frequency axis
cutoff=ceil(N/2); % cut the first half of FFT
Xraw=2*Xraw(1:cutoff)/N;
phase=phase(1:cutoff);
max_nr_harm=7; % in practice the max number of harmonics is determined by the Nyquist frequency
I_n=zeros((max_nr_harm+1)/2,1); % amplitudes of the odd harmonics
delta_n=zeros((max_nr_harm+1)/2,1); % phase angles of the odd harmonics
% extract the harmonic information
for n=1:(max_nr_harm+1)/2
I_n(n)=abs(Xraw(round((2*n-1)*fi/freq_res)+1));
delta_n(n)=phase(round((2*n-1)*fi/freq_res)+1)+0.5*pi; % for the phase angles to be correct 0.5*pi must be added. The reason is unknown.
end
DC=abs(Xraw(1))*cos(angle(Xraw(1)))/2; % Normalization of the DC component need not multiply with 2, so the division by 2. The cosine factor determine the sign of DC offset.
fs2=300; % reconstruct the signal at a higher sampling frequency
t2=0:1/fs2:t(end);
t2=t2';
reconst=zeros(length(t2),2);
reconst(:,1)=sin(2*pi*fi*t2); % the input signal at high sampling frequency
% reconstruct the output signal by summing up the sine components
for n=1:2:max_nr_harm
reconst(:,2)=reconst(:,2)+I_n((n+1)/2)*sin(n*2*pi*fi*t2+delta_n((n+1)/2));
end
reconst(:,2)=reconst(:,2)+DC; % DC shift
plot(x(:,1),x(:,2),'k',reconst(:,1),reconst(:,2),'r')
为了展示，该代码生成一个示范信号。这个示范信号的频率、取样频率和长度都尽可能的任意（频率是个除不尽的小数，取样频率是个质数，长度是12345），也增加了比较明显的噪音。波形重构时，频率是通过对输入信号进行正弦波拟合来获得的。这一步对原信号的噪音还是十分敏感的。如果在以上代码给原信号加上高于千分之一的噪音，整个代码的结果就完全不及格了。这就要由于是后面决定谐波所处的frequency bin那步对信号频率的准确性十分敏感。
有了信号的频率，就可以进行相干取样、FFT……一切都按正常进行。
以上代码本身会给出一个结果比较。我用以上方法对一个真实的LAOS信号进行重构，得到的结果如下(黑的是实验结果，红的是重构曲线）：
实际信号示范
这是一个比较稀有的记录片。斯坦威相比演奏家之潮兴起的历史而言，是一个相对年轻的钢琴厂。对于这些钢琴厂而言，受到大钢琴家的赏识并且在音乐会上使用，是生死攸关的宣传工作。因此这些钢琴厂都很注意搞好与大钢琴家的关系。当然，重要的还是实力，斯坦威的制造工艺具有某种征服性的作用，几乎一下子所有欧洲钢琴家都使用斯坦威了。
什么时候，仪器厂商跟实验家的关系能像这样就好了。
Rubinstein在影片里说，不赞成同一晚上演奏肖邦的两个协奏曲，它们太相似了。也不赞成把Rachmaninoff的两个协奏曲（2和3吧）放一起演。然后开始炫耀他收藏的肖邦手稿、勃拉姆斯私信等。不爱用hobby这词要用passion。
影片是德国人拍摄的。由于二战的关系，犹太人Rubinstein说过不会再去德国演奏。整个影片的拍摄似乎不舍得放过任何一个画面似的，片中的人对Rubinstein也充满了崇敬。在片子中，Rubinstein弹过几个作品的片段，没有弹过完整的，这些片段的拍摄都显得特别认真，旁人马上闭嘴，镜头马上放在最优角度。每当Rubinstein中断演奏，都能强列感到在场者的遗憾。在影片最后，旁述还不无沮丧地说，可惜Rubinstein再也不会在我国公开演奏了！
我们常常说，日本对自己在二战时的罪行反思不足，很多愤青也常常表示要抵制日货。可是，有多少日本人会为中国任何领域的任何人对日本的抵制感到可惜呢？中国人有哪位能像鲁宾斯坦那样为其祖国挣到这样的面子呢？