
f=4000:4000:1100000;f0=2.208*10^6;Kg=.1104
fh=1.104*10^6;
fl=4000;
fh1=f0/2;
fh=25875;
alfa=11.96;np=10;
alfa1=7.09;
m=length(f)
hpf=ones(size(f));
si=ones(size(f));
lpf=ones(size(f));

for i=1:length(f)
hpf(i)=(f(i)^alfa1+fl^alfa1)/(f(i)^alfa1+fh^alfa1);
si(i)=(sin(pi*f(i)/f0)/(pi*f(i)/f0))^2;




lpf(i)=fh1^alfa/(f(i)^alfa+fh1^alfa);


psd(i)=Kg*(2/f0)*si(i)*hpf(i)*lpf(i);
psdnext(i)=psd(i)*.8536*10^(-14)*f(i)^1.5*np^(.6);
psnext(i)=10*log10(1000*psdnext(i));
ps(i)=10*log10(1000*psd(i));
end
plot(f,ps)
%axis([0 f0/2 -150 -85]);grid
hold on
plot(f,psnext)
f1=4000:4000:40000;
hpf1=hpf(1:1,1:10);
gtext('densitatea spectrala de putere (dsp) in aval');
gtext('dsp pentru paradiafonie: aval in amonte');

