function X=mySoundPS(filename,wtype,period)
stacksize(20000000);
L=int(22050*period);
y=read(filename,-1,1)';
N=length(y);
loop=int(N/L);
w=window(wtype,L);
t=linspace(0,period,L+1);
t=t(1:L);
Df=1/period;
f=Df*[-(L/2):(L/2)-1];
fmax=Df*(L/2);
zz=zeros(1,L);
for i=0:loop-1
v=w . *y((i *L+1):(i *L+L));
u=fft(v,-1) ./L;
z=abs(u) .^2;
z=[z(int(L/2)+1:L) z(1:int(L/2))];
zz=zz+z;
power_dB=10*log10(z);
xbasc();
subplot(2,1,1); xgrid(4);
r=[i*period,-1,(i+1)*period,1];plot2d(t+(i*period),y((i*L+1):(i*L+L)),[5],rect=r);
xtitle('Waveform','Time[sec]');
subplot(2,1,2); xgrid(4);
rr=[-fmax,-80,fmax,0];
plot2d(f,power_dB,5,rect=rr);
xtitle('Power[dB]','[Hz]');
printf("time=%1.2f--%1.2f [sec]\n",i*period,(i+1)*period);
playsnd(y((i*L+1):(i*L+L)));
xclick();
end
zz=zz ./loop;
xbasc(); xgrid(4);
power_dB=10*log10(zz);
plot2d(f,power_dB,5,rect=rr);
xtitle('Average Power [dB]','[Hz]');
X=[f' power_dB'];
endfunction
