function handle = wavspectrum( fname, roffset, rduration, q, dolog)

%+++++++++++++++++++++++++++++++++++++
%++ Reading digital sound

N=30;

if nargin < 5
  dolog = 0;
end
if nargin < 4
  q = 0;
  if nargin < 3
    %-- sound duration to consider for the plot (in s)
    rduration = 0.5;
    %-- offset to skip for the leading garbage (in s)
    if nargin < 2 
      roffset = 1;
    end
  end
end

[sounds, fs, bits]=wavread( fname ); 
sound = sounds(:,1);
sound = sound - mean(sound);
slength = length(sound);
disp(['Processing ' fname ]);
disp(['            ' num2str( bits) ' bits']);
disp(['            ' num2str( fs ) ' Hz']);
disp(['            ' num2str( slength / fs) ' s']);

%++
%++ If fs < 44.1kHz, resample to bring the frequency to 44.1kHz
%++
p = round(44100 / fs);
%p = -1;
%%if p > 1
%%  fs = 44100;
%%  sound = resample( sound, p, 1, N);
%%  slength = length(sound);
%%end

%++
%++ Skipping the first second to avoid the impact noise of the instrument
%++
disp(['Skipped: ' num2str( roffset) ' s']);
offset = roffset * fs;
%-- maximum power of 2 fitting in the signal
p2maxl = 2^( fix( log( slength-offset ) / log(2)) );

rlength = min( 2^round( log( rduration * fs) / log(2)), p2maxl );
disp(['Considered duration: ' num2str( rlength/fs) ' s  ( ' num2str( rlength ) ' samples)']);
rsound =sound(offset : offset+rlength);
rslength = length( rsound);

if p > 1
  sscale=1:3*rslength;
else
  sscale=1:rslength;
end
sscale=sscale/rslength * fs;
fsound=abs(fft( rsound ));

%-- plot the aliasing filter impulse response?
if q > 0 
  disp(['The new sampling frequency would have been ' num2str( fs/q) 'Hz']);
  h = antialiasingfilter( q, rslength, N);
end

maxampfsound=max(fsound);
fsound = fsound/maxampfsound;

if dolog > 0
  handle=semilogy(sscale(1:fix(rslength/2)), fsound(1:fix(rslength/2)), 'b');
  hold on
  if p > 1
    handle=semilogy(sscale,[ fsound' fsound' fsound' ]','r');
  end
  handle=semilogy(sscale(1:fix(rslength/2)), fsound(1:fix(rslength/2)), 'b');
else
  handle=plot(sscale(1:fix(rslength/2)), fsound(1:fix(rslength/2)),'b');
  hold on
  if p > 1
    handle=plot(sscale, [ fsound' fsound' fsound' ]','r');
  end
  handle=plot(sscale(1:fix(rslength/2)), fsound(1:fix(rslength/2)),'b');
end
hold off
if q > 0
  hold on
  if dolog > 0
    semilogy( [fs/q fs/q], [1e-8 1.05],'k');
    semilogy( [fs/q/2 fs/q/2], [1e-8 1.05],'g');
    semilogy( sscale, abs(fft(h)) ,'r');
  else
    plot( [fs/q fs/q], [0 1.05],'k');
    plot( [fs/q/2 fs/q/2], [0 1.05],'g');
    plot( sscale, abs(fft( h)) ,'r');
  end
  hold off
end
xlabel('\omega (Hz)');
ylabel('|f(\omega)|');
if dolog > 0 
  axis([0 5000 1e-4 1.05]);
else
  axis([0 5000 0 1.05]);
end