Thread: POSSIBLE error in trans.c

Dear Mrs. / Miss,

First, please allow me to introduce myself. I'm Youshan Liu, a student at Institute of Geology and Geophysics, Chinese Academy of Sciences.

Recently, I found several problems in SAC.

1. On lines 84-86 in trans.c


/* Set the real part of the Frequecy Response at the Nyquist Frequency to
the Amplitude of the Nyquist Freuquency */


sre[nfft - 1] =
sqrt(sre[nfft - 1] * sre[nfft - 1] + sim[nfft - 1] * sim[nfft - 1]);
sim[nfft - 1] = 0.0;




I think there are two errors. 1) sre[nfft-1] is not the Nyquist Frequency point, it should be sre[nfft/2]. 2) Why a sqrt is done. When it is negative. Sure, the Fourier cofficient at the Nyquist frequency should be real instead of complex, while it can be negative float number.

As the above presented, it must be positive.

2. In lines 249-259 in transfer.c


for (i = 1; i < nfreq; i++) {
denr = (powi(sre[i], 2) + powi(sim[i], 2));
if (denr <= FLT_MIN) {
sre[i] = 0.0e0;
sim[i] = 0.0e0;
} else {
denr = 1.0e0 / denr;
sre[i] = sre[i] * denr;
sim[i] = -sim[i] * denr;
}
}





I think that it is not a good idea to use FLT_MIN=1.175494350822287507968736537222245678e-38 to stablize deconvolution. It should be data adaptive, such as water_level * max(abs(sre*sre + sim*sim)).

In addition, it may be implemented as follows.

Assuming x is the real waveform, h is the response of seismometer, while y is the observed data. They have the following relationship, x * h = y.

Thus, x = y / h. A stable version of it is x = conj(h) * y / max(conj(h)*h, wtr), where wtr = water_level * max(conj(h)*h).




I'm looking forward your reply!




Thanks a lot!




Best Wishes,

Youshan Liu