function []=plot_spectrum(data,t1,t2)
f=1/(data(2,1)-data(1,1));
to=data(1,1);
n1=floor((t1-to)*f)+1;
n2=floor((t2-to)*f)+1;
data=data(n1:n2,:);
t=data(:,1);
% Plot the time series of the wave form
% figure()
% plot(t,y,'r-');
% xlabel('Time (s)');
% ylabel('Signal amplitude (Pa)');

% Calculate the spectrum using Welch's method

zeros = 2048*8;                   % number of zero padded cells
nfft = round(length(data(:,1))+zeros-1); % FFT length with zero padding
nfft = 2^nextpow2(nfft);          % FFT length with zero padding
lwin = length(data(:,1));                % Window size for Welch's method, use the length of the original signal
Hzn = f;                     % Sampling frequency

% psd_y=zeros(nfft,6);
% F=zeros(nfft,6);

for k=2:5
y=data(:,k);
[psd_y(:,k-1),F(:,k-1)] = pwelch(detrend(y),lwin,[],nfft,Hzn);
end

% Plot the spectrum

figure()
% plot(F(:,1),psd_y(:,1),F(:,2),psd_y(:,2),F(:,3),psd_y(:,3),F(:,4),psd_y(:,4),F(:,5),psd_y(:,5),F(:,6),psd_y(:,6));
plot(F(:,1),psd_y(:,1),F(:,2),psd_y(:,2),F(:,3),psd_y(:,3),F(:,4),psd_y(:,4));
xlim([0 30])
title(['Power Spectrum for all sensors, ',num2str(t1),' to ',num2str(t2),'s']);
xlabel('Frequency (Hz)');
ylabel('Spectral Density (Pa^2/s)');
 legend('1','2','3','4','5','6')
end