function freqs = image_freq(f,fs,fmin,fmax)
% IMAGE_FREQ produces the imaged frequencies that correspond to
% the frequency and sampling frequency given
%
% FREQS = IMAGE_FREQ(F,FS,FMIN,FMAX) 
%
% returns the array of image frequencies FREQS that lie between FMIN and FMAX
% using k*FS +/- F

if(f < 0 | fs < 0)
   error('Negative frequencies not allowed for signal frequency or sampling frequency.');
end
kmin=ceil((fmin-f)/fs);          % find lowest necessary index
kmax=floor((fmax+f)/fs);         % find highest necessary index
i=1;
for k=kmin:kmax                  % generate all possible frequencies
      freqs(i) = k*fs-f;
      freqs(i+1) = k*fs+f;     
      i= i+2;
end
freqs = sort(freqs);   % sort
i=1;
while(freqs(i) < fmin)      % weed out frequencies outside the specified range
   i = i+1;
end
n1=i;
i=length(freqs);
while(freqs(i) > fmax)
   i = i-1;
end
n2=i; 
freqs=freqs(n1:n2);