5
$\begingroup$

Equation (4) of a recent paper (An improved peak-frequency-shift method for Q estimation) shows how to generate a wavelet using a zero-phase Ricker wavelet.

$|B_1(f)| = |B_0(f)|\exp(-\pi t f / Q)$

In the above, $|B_1(f)|$ is the amplitude spectra of the wavelet that I want to calculate, $|B_0(f)|$ is the amplitude spectra of the reference (Ricker wavelet), $\pi$ is the ubiquitous constant, $t$ is the time separation between wavelets, $f$ is the frequency, and $Q$ is a number (the Quality factor).

I know how to calculate the zero-phase Ricker wavelet in the time domain as $B_0(t)$, and I know how to calculate $|B_0(f)|$ as the square root of the sum of the real and imaginary parts squared of the frequency-domain reference Ricker wavelet. The real and imaginary parts are found using the FFT.

However, now that I've calculated $|B_1(f)|$, how do I transform $|B_1(f)|$ back into the time domain? I know $|B_1(f)|$, but can I find the real and imaginary parts of $|B_1(f)|$ so that I can use the IFFT to compute $B_1(t)$ in the time domain, or is some other information missing and required?

Is the "without phase changes" words in the recent paper (http://www.seg.org/documents/10161/997244/1218.pdf) a clue that might help with the numerical implementation?

So how did the authors generate the wavelets in Figure 1 of the paper? Is something other than Equation (4) being used here?

Perhaps there is nothing more happening here than an issue with notation. Looking at another paper (Estimation of quality factors from CMP records), the frequency spectrum of a Ricker wavelet is computed using Equation (1):

$B(f) = \frac{2}{{\sqrt \pi }}\frac{{{f^2}}}{{f_m^2}}{e^{\frac{{ - {f^2}}}{{f_m^2}}}}$

In the above, $f_m$ is the dominant frequency. However, there doesn't seem to be any complex part of this expression. Does this imply that the complex part is comprised of all zeros? To numerically evaluate this equation, how do I compute the $f$ variable? For example, using Matlab, $f$ can be computed using the colon operator:

f = fmin:deltaf:fmax

But what is fmin and fmax? Now $B(f)$ can be computed and taken back into the time domain using the IDFT (or the IFFT), but what if I want a real time-domain discrete signal for the $B(t)$ function?

Equation (5) of the same paper gives an expression that I think can be used to calculate the Q factor wavelet (note that this looks the same as the equation at the very top of this post):

$B(f,t) = B(f){e^{\frac{{ - \pi ft}}{Q}}}$

However, what is meant here by the "amplitude spectrum"? How do I get the real and complex parts of this spectrum?

A time domain version equation of the Ricker wavelet can be found here: Wikipedia link, but then what is the time domain expression used to calculate the Q factor wavelet? As an aside, it would have been easier to have the time domain version of:

$B(f,t) = B(f){e^{\frac{{ - \pi ft}}{Q}}}$

But what is the "amplitude spectrum" in this context, and how can I compute the wavelet in the time domain?

$\endgroup$
3
$\begingroup$

The complex result of a DFT of a signal codifies amplitude and phase information. When you take a magnitude of the resulting complex vector, you automatically throw away any phase information.

In so many words, the DFT result is informing you of what amplitudes and phases of all its frequency components are required to reconstruct your signal. If you strip away phase information you are lining up all your complex sinusoids and then adding them up coherently. If you lose phase how can you reconstruct your signal?

That being said, what you can get in this case is the autocorrelation function of $b_1(t)$, by taking the IDFT of $|B_1(f)|^2$, but you need phase information to completely reconstruct the original signal.

$\endgroup$
  • 1
    $\begingroup$ The autocorrelation function of a signal actually transforms to its power spectral density, so you would need to take the IDFT of $|B_1(f)|^2$. $\endgroup$ – Jason R May 17 '12 at 1:52
  • 1
    $\begingroup$ @JasonR Oops - thanks for that - have corrected. $\endgroup$ – Spacey May 17 '12 at 5:12
  • $\begingroup$ @Mohammad: Thank you; this does indeed make sense. If this is the case, then how did the authors of the paper "generate a series of target wavelets" using Equation (4)? There is something else being done other than evaluation of Equation (4)? $\endgroup$ – Nicholas Kinar May 17 '12 at 13:33
  • $\begingroup$ @NicholasKinar I am not sure, I will have to read the paper and get back to you on it. $\endgroup$ – Spacey May 17 '12 at 15:20
  • $\begingroup$ @Mohammad: OK, thank you. I've also added some additional information in my question above. $\endgroup$ – Nicholas Kinar May 17 '12 at 15:25
1
$\begingroup$

The issue with calculating the wavelet appears to be nothing more than an issue with notation and terminology. Starting with the second paper mentioned in my original post (Estimation of quality factors from CMP records), the magnitude (real part) of the Ricker wavelet in the frequency domain is given by:

$B(f) = \frac{2}{{\sqrt \pi }}\frac{{{f^2}}}{{f_m^2}}{e^{\frac{{ - {f^2}}}{{f_m^2}}}}$

In the above, $\pi$ is the ubiquitous constant (3.1415...), whereas $f$ is a frequency vector ranging between two frequencies (0 Hz and 400 Hz for example), and $f_m$ is the center frequency of the wavelet. The phase (imaginary part) in the frequency domain is a vector comprised of zeros.

The Ricker wavelet at a time difference $t$ and a known $Q$ is calculated using:

$B(f,t) = B(f){e^{\frac{{ - \pi ft}}{Q}}}$

These equations are used to generate a filter kernel in the frequency domain. After taking the IDFT (or IFFT), the discrete time domain signal needs to be shifted (see Chapter 17 of Smith), but I don't think that I need to truncate and window the kernel, since this will change the frequency response.

I think that vector Bf2 in the code below needs to be symmetric to preserve the symmetry of the FFT, and I've calculated deltat from fmax.

A paper that discusses Q estimation can be found here (Estimation of Q using crosscorrelation).

The following Matlab code snippet demonstrates what I've done:

% inputs
fm = 50;            % dominant frequency of wavelet (Hz)
deltaf = 0.1;       % frequency step (Hz)
fmax = 400;         % maximum frequency (Hz)
tdiff = 0.5;        % time propagation distance of wavlet (s)
Q = 40;             % Q value of wavelet

deltat = 1 / (2 * fmax);
f = 0:deltaf:fmax;
term = (f.^2) ./ (fm^2);
term1 = pi .* tdiff .* f ./ Q;
Bf = (2 / sqrt(pi)) .* (term) .* exp(-term) .* exp(-term1);

figure; plot(f, Bf);
title('Frequency spectrum of Ricker wavelet')
xlabel('Frequency (Hz)')
ylabel('Amplitude')

% We could convert the magnitude and phase (all zero) into rectangular
% coordinates from polar coordinates
% This doesn't really do anything in this case, since the phase
% is just a zero vector
% mag = real(Bf);
% phase = imag(Bf);
% Rp = mag .*  cos(phase);
% Ip = mag .* sin(phase);

%  prepare the array for taking back into the time domain
%  this should be symmetric around the N/2 + 1 point
Bf2 = [Bf fliplr(Bf(1:end-1))];

% take back into time domain 
Bt =  ifft(Bf2, 'symmetric');

% rotate
N = length(Bt);
Btr = zeros(N,1);
for i = 0:N-1
   index = floor(i + N/2);
   if (index > N-1)
    index = index - N; 
   end
   Btr(index+1) =  Bt(i+1); 
end

% determine the time axis
t = deltat:deltat:deltat*(N/2);
tt = [-fliplr(t) 0 t]; 

figure; plot(tt, Btr);
title('Time spectrum of Ricker wavelet')
xlabel('Time (s))')
ylabel('Amplitude')
xlim([-0.08 0.08])
$\endgroup$
  • $\begingroup$ Does this appear to be reasonable (particularly the symmetry of the Bf2 vector)? $\endgroup$ – Nicholas Kinar May 18 '12 at 2:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.