function [rss,rssfield]=rssound(filename,spec,temp,f0,f1,dur,fs);
%RSSOUND generates a random spectrogram sound
%
%   RSSOUND(filename,spec,temp);
%   generates a random spectrogram sound with temp temporal components
%   and spec spectral components. The sound is written into the file 
%   specified by filename.
%   
%   RSSOUND(filename,spec,temp,f0,f1,dur,fs);
%   generates a random spectrogram sound with temp temporal components
%   in dur seconds (if dur = 1, temp is the modulation frequency in Hz)
%   and spec spectral components in the frequency band between
%   f0 and f1 (in Hz). The sound is written into file filename with
%   sampling frequency fs.
%
%   [rss,rfield] = RSSOUND(filename,...);
%   returns the samples of the random spectrogram sound in rss
%   and the initial random field used to generate the sound in rfield.
%   
%   If f0,f1,dur or fs are ommited, the following default values are used:
%   f0 = 250Hz, f1 = 8kHz, dur = 1sec, fs = 44.1kHz.
%
%   All amplitude changes within the sound are gated with 10 ms linear ramps, 
%   the sound on- and offset are gated with 10ms cosine ramps.
%
%   Example:
%           [rss,rssfield]=rssound('rss3x5',3,5,500,4000,1,44100);
%           imagesc(rssfield); % show the initial random field
%           sound(rss,44100); % play the random frequency sound
%
%	This script is published as supplementary material to the article by
%   Schönwiesner M, Rübsamen R and Cramon DY. "Hemispheric Asymmetry for
%   Spectral and Temporal Processing in the Human Antero-lateral Auditory
%   Belt Cortex" submitted to the European Journal of Neuroscience.
%
%
%   (c) 2004 Marc Schönwiesner
%   $Revision: 1.0 $  $Date: 2004/05/20 15:37:25$

% check arguments
if ( nargin < 3 )
  help rssound
  return
end

% standard values:
dur_std = 1;        % standard value for len (1 second)
f0_std  = 250;      % standard value for f0 (250 Hz)
f1_std  = 8000;     % standard value for f1 (8 kHz)
fs_std  = 44100;    % standard value for fs (44.1 kHz)
numSinPerErb = 50;

% check dur
if ~(exist('dur','var'))
  dur = dur_std;
elseif isempty(dur)
  dur = len_std;
end

% check f0
if ~(exist('f0','var'))
  f0 = f0_std;
elseif isempty(f0)
  f0 = f0_std;
end

% check f1
if ~(exist('f1','var'))
  f1 = f1_std;
elseif isempty(f1)
  f1 = f1_std;
elseif (f1 < f0)
  error('f1 < f0')
end

% check fs
if ~(exist('fs','var'))
  fs = fs_std;
elseif isempty(fs)
  fs = fs_std;
elseif (fs < 2*f1)
  error('fs too low for specified upper frequency band')
end

% set variables
spec = ceil(abs(spec));
temp = ceil(abs(temp));
numErb = (Hz2Erb(f1) - Hz2Erb(f0));   % number of ERBs between f0 and f1
numSin = ceil(numErb * numSinPerErb);
numSinPerSpec = numSin / spec;
numSamples = fs * dur;

% generate the initial random field
% please constrain the random fields according to your requirements by
% adding code here
rssfield = rand(spec,temp);

% please constrain the random fields according to your requirements by
% adding code here (for instance to keep the stimulus energy within certain bounds)
rssfield(rssfield < 0.5) = 0;
energy = sum(sum(rssfield)')/(spec*temp);

% construct linear ramps of 10 ms around amplitude changes
disp('ramping of amplitude changes...');
rsslice = zeros(spec,numSamples);
mult = round(dur*100/temp);
xi = linspace(1,(mult*temp),numSamples);
for I = 1:spec
    tmp = repmat(rssfield(I,:),mult,1);
    tmp = tmp(:);
    tmp = interp1(1:(mult*temp),tmp(:),xi,'linear');
    rsslice(I,:) = tmp;
end

disp('computing ERB Frequencies...');
cfs=computeErbFs(f0,f1,numSin);

disp('generating sound...');
ivec = linspace(spec+1,1,numSin);
ivec(1) = spec; ivec = floor(ivec);
rss=zeros(1,numSamples);
for I = numSin:-1:1
    intslice = rsslice(floor(ivec(I)),:);
    rss = rss + (msPureTone(cfs(I),dur*1000,fs).*rsslice(floor(ivec(I)),:));
    fprintf(1,'%d|',I);
end

fprintf(1,'\nDone.');
% gating the sound
gate = fs / 100;
rss(1:gate) = rss(1:gate).*sin((0:gate-1)*pi/(2*(gate-1))).^2;
rss(end-gate+1:end) = rss(end-gate+1:end).*cos((0:gate-1)*pi/(2*(gate-1))).^2;

% output
rss=rss*0.85/max(rss);
wavwrite(rss,fs,filename);


% end

% subfunctions
function y = computeErbFs(f0,f1,num)
y=Erb2Hz(linspace(Hz2Erb(f0),Hz2Erb(f1),num));

function y=Erb2Hz(x)
y=(10.^(x/21.4)-1)/4.37e-3;

function y=Hz2Erb(x)
y=(21.4*log10(4.37e-3*x+1));

function x=msPureTone(freq,dur,fs,phase,level)
sdur=round(fs*dur/1000);
level=80;ref=80;
phase = rand*2*pi*freq;
amp=10.^((level-ref)./20);
x=amp.*sin(2*pi*(freq/fs)*((1:sdur)+phase));