function han = sigspec(sig,NFFT,fs,flag)
% sig  = input signals in the columns
% NFFT = number of points in the fast Fourier transform
% fs   = sample rate
% flag = 'lin' for linear magnitude plot, 'log' for log magnitude plot

if(nargin < 2)
NFFT = 2^(ceil(log2(size(sig,1)))+1);
end
if(nargin < 3)
fs = 1;
end
if(nargin < 4)
flag = 'lin';
end

f=[0:NFFT-1]/NFFT; % Normalized frequency vector
f = f - 0.5; % Shifted normalized frequency vector
f = f*fs; % Properly scaled frequency vector

SIG = fftshift(fft(sig,NFFT)); % Compute discrete Fourier transform

if(strcmp(flag,'lin'))
han = plot(f,abs(SIG)); % Plot magnitude (linear scale)
else(strcmp(flag,'log'))
han = plot(f,20*log10(abs(SIG))) % Plot magnitude (dB scale)
end
set(gca,'Xlim',[-fs/2 fs/2]);