A nice combination of NumPy and ObsPy enables us to quickly plot FFTs from the event SAC files. Finding FFT is the easy part, thanks to the 'numpy.fft' method. The main roadblock for us, as it turned out, was the x-axis scale part.
To plot a FFT, matplotlib.pyplot needs to be provided with two arrays, sets of y points and corresponding x points. NumPy.fft does the former for us. Latter can be created through another method of Numpy: 'fft.fft.fftfreq (n=number of points,d=spacing in sample)'.
We were able to locate this function quickly, however we committed mistakes in implementing it, all because we were giving it the wrong arguments. The documentation declares 'd' as spacing, which we thought was the spacing between two consecutive data points along the x-Axis in the frequency domain.
// d= sampling_freq/total Number of points
This turned out to be wrong, and gives incorrect domain to the the DFT, and the scale showed ticks above Nyquist frequency. So we looked at the source code, and found that 'd' represented distance between two sample in the time domain.
// d= 1/sampling_freq
This gave us the x-scale that was expected, i.e. within the range of Nyquist frequency. Below is one such plot showing the correct scale, for all the three channels.
To plot a FFT, matplotlib.pyplot needs to be provided with two arrays, sets of y points and corresponding x points. NumPy.fft does the former for us. Latter can be created through another method of Numpy: 'fft.fft.fftfreq (n=number of points,d=spacing in sample)'.
We were able to locate this function quickly, however we committed mistakes in implementing it, all because we were giving it the wrong arguments. The documentation declares 'd' as spacing, which we thought was the spacing between two consecutive data points along the x-Axis in the frequency domain.
// d= sampling_freq/total Number of points
This turned out to be wrong, and gives incorrect domain to the the DFT, and the scale showed ticks above Nyquist frequency. So we looked at the source code, and found that 'd' represented distance between two sample in the time domain.
// d= 1/sampling_freq
This gave us the x-scale that was expected, i.e. within the range of Nyquist frequency. Below is one such plot showing the correct scale, for all the three channels.
 |
| Stacked FFTs, with sampling rate of 50Hz. Notice the max freq on X axis is half sampling rate, or the Nyquist frequency |
No comments:
Post a Comment