利用arma、ar、ma模型，以及周期图等进行系统参数估计.m

利用ARMA、AR、MA模型

clear

tic
f1 = 1776;
f2 = 1000;
fs = 44100;
n = 1:4410;
x = sin(2*pi*f1*n/fs) + sin(2*pi*f2*n/fs);

IMF = emd(x, 'MAXMODES', 3);
        
lay = abs(fft(IMF(1,:))); 
[ pk pos ] = max(lay(1:end/2));
lay_f = (pos-1)*10;   
fm_est1 = xray_cal_f(IMF(1,:), fs);
buf1 = sprintf('IMF(1,:)主要频率成分是%fHz, 估计得到的MASK信号频率是%fHz',lay_f, fm_est1);
disp(buf1);

lay = abs(fft(IMF(2,:))); 
[ pk pos ] = max(lay(1:end/2));
lay_f = (pos-1)*10;   
fm_est2 = xray_cal_f(IMF(2,:), fs);
buf2 = sprintf('IMF(2,:)主要频率成分是%fHz, 估计得到的MASK信号频率是%fHz',lay_f, fm_est2);
disp(buf2)

xm1 = sin(2*pi*fm_est1*n/fs);
IMF = emd(x, 'MAXMODES', 3, 'MASK', xm1);
lay = abs(fft(IMF(1,:))); 
[ pk pos ] = max(lay(1:end/2));
lay_f1 = (pos-1)*10;   
lay = abs(fft(IMF(2,:))); 
[ pk pos ] = max(lay(1:end/2));
lay_f2 = (pos-1)*10;   
buf1 = sprintf('采用%fHz的MASK信号, 提取出来的IMF(1,:)是%fHz信号, IMF(2,:)是%fHz信号',fm_est1,lay_f1,lay_f2);
disp(buf1);
figure(1);hold on;
plot(IMF(1,:),'r');
plot(IMF(2,:),'b');
xm2 = sin(2*pi*fm_est2*n/fs);
IMF = emd(x, 'MAXMODES', 3, 'MASK', xm2);
lay = abs(fft(IMF(1,:))); 
[ pk pos ] = max(lay(1:end/2));
lay_f1 = (pos-1)*10;   
lay = abs(fft(IMF(2,:))); 
[ pk pos ] = max(lay(1:end/2));
lay_f2 = (pos-1)*10;   
buf2 = sprintf('采用%fHz的MASK信号, 提取出来的IMF(1,:)是%fHz信号, IMF(2,:)是%fHz信号',fm_est2,lay_f1,lay_f2);
disp(buf2);
figure(2);hold on;
plot(IMF(1,:),'r');
plot(IMF(2,:),'b');

disp('结论：MASK信号的频率选择，对IMF信号的提取有很大影响！');

toc