求助从一个叠加了扰动量的正弦信号中提取出扰动量的MATLAB程序

血色黄昏

    部分程序如下，我看的不是很懂。具体情况是一个电力系统模型稳定运行后，20s在发电机侧添加一个振荡扰动，时长6s，扰动消失后系统重新稳定。在workspace中保存了线路中一个节点采样得到的电压值，如u1，采样点为60000个，现在要提取u1中在添加扰动后叠加的扰动量，即通过u1电压值分离得到扰动量的波形。右图中上部分为u1，稳定时为一个正弦信号。有没有哪位高手能帮忙看看这个问题，优化下程序，万分感激，先谢谢了。
















plot (u1(1:60000,1), 'DisplayName', 'u1(1:60000,1)', 'YDataSource', 'u1(1:60000,1)'); figure(gcf)
a1=hilbert(u1)
hu1=abs(a1)
plot (hu1(1:60000,1), 'DisplayName', 'hu1(1:60000,1)', 'YDataSource', 'hu1(1:60000,1)'); figure(gcf)
a2=hilbert(u2)
hu2=abs(a2)
plot (hu2(1:60000,1), 'DisplayName', 'hu2(1:60000,1)', 'YDataSource', 'hu2(1:60000,1)'); figure(gcf)
a3=hilbert(u3)
hu3=abs(a3)
plot (hu3(1:60000,1), 'DisplayName', 'hu3(1:60000,1)', 'YDataSource', 'hu3(1:60000,1)'); figure(gcf)
plot (hu2(1:60000,1), 'DisplayName', 'hu2(1:60000,1)', 'YDataSource', 'hu2(1:60000,1)'); figure(gcf)
plot (hu1(1:60000,1), 'DisplayName', 'hu1(1:60000,1)', 'YDataSource', 'hu1(1:60000,1)'); figure(gcf)
hju1=hu1(30000:50000)
mm1=max(hju1)
mi1=min(hju1)
hjmu1=hju1-(mm1-mi1)/2-mi1
plot (hjmu1(1:20001,1), 'DisplayName', 'hjmu1(1:20001,1)', 'YDataSource', 'hjmu1(1:20001,1)'); figure(gcf)
hju3=hu3(30000:50000)
mm3=max(hju3)
mi3=min(hju3)
hjmu3=hju3-(mm3-mi3)/2-mi3
plot (hjmu3(1:20001,1), 'DisplayName', 'hjmu3(1:20001,1)', 'YDataSource', 'hjmu3(1:20001,1)'); figure(gcf)
plot (hjmu1(1:20001,1), 'DisplayName', 'hjmu1(1:20001,1)', 'YDataSource', 'hjmu1(1:20001,1)'); figure(gcf)
a=xcorr(hjmu3,hjmu1)
[m13,t13]=max(a)
plot (hjmu1(1:20001,1), 'DisplayName', 'hjmu1(1:20001,1)', 'YDataSource', 'hjmu1(1:20001,1)'); figure(gcf)
hjmju1=hjmu1(2500:12500)
plot (hjmju1(1:10001,1), 'DisplayName', 'hjmju1(1:10001,1)', 'YDataSource', 'hjmju1(1:10001,1)'); figure(gcf)
hjmju3=hjmu3(2500:12500)
plot (hjmju3(1:10001,1), 'DisplayName', 'hjmju3(1:10001,1)', 'YDataSource', 'hjmju3(1:10001,1)'); figure(gcf)
plot (hjmju1(1:10001,1), 'DisplayName', 'hjmju1(1:10001,1)', 'YDataSource', 'hjmju1(1:10001,1)'); figure(gcf)
a=xcorr(hjmju1,hjmju3)
[m13,t13]=max(a)
plot (a(1:20001,1), 'DisplayName', 'a(1:20001,1)', 'YDataSource', 'a(1:20001,1)'); figure(gcf)
r13=a
plot (r13(1:20001,1), 'DisplayName', 'r13(1:20001,1)', 'YDataSource', 'r13(1:20001,1)'); figure(gcf)
rr13=r13/100000
plot (rr13(1:20001,1), 'DisplayName', 'rr13(1:20001,1)', 'YDataSource', 'rr13(1:20001,1)'); figure(gcf)

yjhou

没看懂，不过我看思路好像是有傅立叶分析。