Find the ratio or difference spectra as time progresses

Hi, I have a number of spectra that has been collected over a certain number of hours. How do I compute the difference spectra or ratio spectra over that time period? So say spectrum 1 is rationed with spectrum 2, spectrum 3 rationed with spectrum 2 etc. can same be done with difference spectrum as well? Any algorithm to do this? Thank you.

Answer 1

[Star Strider]

If the spectra are vectors and have a relatively flat baseline (I am not certain what your spectra are, since they can be mass spectra, NMR spectra, Fourier taransforms of time-domain signals, among others),and if they all share the same independent variables (for example a specific range of frequencies), taking their differences and ratios would be relattvely straightforward.

x = linspace(0, 1E+3, 1E+4);
p = randi(1E+3, 5, 1);
y = exp(-(x-p).^2/100)./rand(size(p))+0.1;

figure
tiledlayout(size(y,1),1)
for k = 1:size(y,1)
    nexttile
    plot(x, y(k,:))
    grid
    xlabel('X')
    ylabel('Y')
    title("Spectrum "+k)
end

dy =  y(1,:);
for k = 1:size(y,1)-1
    dy = dy - y(k+1,:); 
end

figure
plot(x, dy)
grid
xlabel('X')
ylabel('Y')
title('Spectra Differences')
axis('padded')

qy =  y(1,:);
for k = 1:size(y,1)-1
    qy = qy ./ y(k+1,:); 
end

figure
plot(x, qy)
grid
xlabel('X')
ylabel('Y')
title('Spectra Quotient')
axis('padded')

It all depends on what your spectra are, and what you want to do with them.

.

Answer 2

[William Rose]

There is a large literature on transfer funciton estimation by taking the ratio of two complex spectra. This usually involves the spectra of two simultaneous signals (input and output). The ratio of output spectrum to input spectrum, better known as the transfer function, describes the properties of the (linear, time invariant) system that transforms the input into the output.

A major result of the analysis of diferent way of estimating the transfer function is that the estimate T(f)=Y(f)/X(f) is NOT a good way to estimate the transfer function. The reasons are somewhat complicated. It turns out that a much better estimate is T(f)=Pxy(f)/Pxx(f), wher Pxx is thre power spectrum of x and Pxy is the cross power spectrum between x and y.

Depending on your situation and goals, you might want to compute Pxy and Pxx for your signals, and take the ratio of those.

Example with signals x and y, and the ratio of their complex spectra: T(f)=Y(f)/X(f). In this example, I simply take the ratio of the spectra for x(t) and y(t). I do not compute the power spectrum or cross power spectrum.

N=1000;             % number of samples
fs=1;               % sampling rate (Hz)
frange=[0.1,0.3];   % bandpass frequency range (Hz)
[b,a]=butter(4,frange/(fs/2));  % bandpass filter
x=randn(1,N);       % x=white noise
y=filter(b,a,x);    % y=bandpass-filtered version of x
% Compute spectrum estimates for x, y
X=fft(x);           % X is complex
X=X(1:N/2+1);       % keep only the lower half of the spectrum 
Y=fft(y);           % Y is complex
Y=Y(1:N/2+1);       % keep lower half of spectrum
T1=Y./X;            % T1=ratio Y/X
% plot the ratio, versus frequency
f=(0:N/2)*fs/N;  % frequency vector (Hz)
figure
subplot(211), plot(f,abs(T1),'-b.'); grid on
title('|T1|=|Y(f)/X(f)|'); ylabel('Amplitude Ratio');
subplot(212), plot(f,angle(T1)*180/pi,'-b.'); grid on
title('Phase(T1)'); xlabel('Frequency (Hz)'); ylabel('Phase (degrees)');

Notice that the amplitude ratio of the spectra (top panel) is approximately 1, from f=0.1 to f=0.3 Hz. This is the pass band of the digital filter. Outside the pass band, the ratio is approximately 0. This is what we expect, since y is a band-passed version of x. The estimates for the ampltude ratio and phase are noisy. The estimates would be less noisy if we computed them with the power spectrum and cross power spectrum, as discussed above.