写一个对信号进行局部均值分解的matlab代码，要求代码可以绘制出原始信号时域波形和每个PF分量的时域波形和包络谱

function [envelope, PF] = LMD(signal, numPF)
% signal: 输入信号
% numPF: 分解出的分量个数

% 设置初始值
h = signal;
envelope = zeros(size(signal));
PF = cell(numPF, 1);

% 局部均值分解
for i = 1:numPF
% 计算局部均值
meanh = movmean(abs(h), 5);
% 计算局部极值点
maxh = islocalmax(abs(h));
minh = islocalmin(abs(h));
% 计算扩展极值点
maxh = max(maxh, meanh > abs(h));
minh = max(minh, meanh > abs(h));
% 计算上下包络线
up = meanh + (h - meanh) .* (1 - maxh);
down = meanh - (meanh - h) .* (1 - minh);
% 计算瞬时频率
instfreq = diff(unwrap(angle(h)));
instfreq = [instfreq(1); instfreq];
% 计算瞬时周期
instperiod = 1 ./ instfreq;
instperiod(isnan(instperiod)) = 0;
% 计算PF分量
PF{i} = h .* instperiod;
% 更新残差
h = h - PF{i};
% 更新包络线
envelope = envelope + up + down;
end

% 绘制时域波形
figure;
subplot(numPF+1,1,1);
plot(signal);
title('原始信号');
ylabel('幅值');
xlabel('时间');
grid on;

for i = 1:numPF
% 绘制PF分量时域波形
subplot(numPF+1,1,i+1);
plot(real(PF{i}));
title(['PF', num2str(i), '分量']);
ylabel('幅值');
xlabel('时间');
grid on;
% 绘制PF分量包络谱
figure;
plot(abs(fft(PF{i})));
title(['PF', num2str(i), '分量包络谱']);
ylabel('幅值');
xlabel('频率');
grid on;
end

% 绘制原始信号包络线
figure;
plot(envelope);
title('原始信号包络线');
ylabel('幅值');
xlabel('时间');
grid on;

end