谁会matlab编程显示功率谱密度呢，跪求...

t=0:0.0001:0.1; %时间间隔为0.0001，说明采样频率为10000Hz
x=square(2*pi*1000*t); %产生基频为1000Hz的方波信号
n=randn(size(t)); %白噪声
f=x+n; %在信号中加入白噪声
figure(1);
subplot(2,1,1);
plot(f); %画出原始信号的波形图
ylabel('幅值(V)');
xlabel('时间(s)');
title('原始信号');
y=fft(f,1000); %对原始信号进行离散傅里叶变换，参加DFT采样点的个数为1000
subplot(2,1,2);
m=abs(y);
f1=(0:length(y)/2-1)'*10000/length(y);%计算变换后不同点对应的幅值
plot(f1,m(1:length(y)/2));
ylabel('幅值的模');
xlabel('时间(s)');
title('原始信号傅里叶变换');
%用周期图法估计功率谱密度
p=y.*conj(y)/1000; %计算功率谱密度
ff=10000*(0:499)/1000; %计算变换后不同点对应的频率值
figure(2);
plot(ff,p(1:500));
ylabel('幅值');
xlabel('频率(Hz)');
title('功率谱密度（周期图法）');