diff --git a/huddle-test-geophones/figs/psd_velocity.png b/huddle-test-geophones/figs/psd_velocity.png index b27c091..4e862f6 100644 Binary files a/huddle-test-geophones/figs/psd_velocity.png and b/huddle-test-geophones/figs/psd_velocity.png differ diff --git a/huddle-test-geophones/figs/tf_geophones.png b/huddle-test-geophones/figs/tf_geophones.png new file mode 100644 index 0000000..eddee9d Binary files /dev/null and b/huddle-test-geophones/figs/tf_geophones.png differ diff --git a/huddle-test-geophones/index.html b/huddle-test-geophones/index.html index 444ad26..814c589 100644 Binary files a/huddle-test-geophones/index.html and b/huddle-test-geophones/index.html differ diff --git a/huddle-test-geophones/index.org b/huddle-test-geophones/index.org index 90dbf4d..78e2066 100644 --- a/huddle-test-geophones/index.org +++ b/huddle-test-geophones/index.org @@ -96,8 +96,8 @@ The voltage amplifiers include a low pass filter with a cut-off frequency at 1kH ** Compute PSD #+begin_src matlab :results none - [pxx1, f1] = pwelch(x1, hanning(ceil(length(t)/100)), 0, [], 1/dt); - [pxx2, f2] = pwelch(x2, hanning(ceil(length(t)/100)), 0, [], 1/dt); + [pxx1, f1] = pwelch(x1, hanning(ceil(1/dt)), 0, [], 1/dt); + [pxx2, f2] = pwelch(x2, hanning(ceil(1/dt)), 0, [], 1/dt); #+end_src ** Take into account sensibility of Geophone @@ -147,6 +147,7 @@ The cut-off frequency is set at 1kHz. set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]') + xlim([2, 500]); #+end_src #+NAME: fig:psd_velocity @@ -159,3 +160,37 @@ The cut-off frequency is set at 1kHz. #+CAPTION: Spectral density of the velocity #+RESULTS: fig:psd_velocity [[file:figs/psd_velocity.png]] +** Transfer function between the two geophones +#+begin_src matlab :results none + [T12, f12] = tfestimate(x1, x2, hanning(1/dt), 0, [], 1/dt); +#+end_src + +#+begin_src matlab :results none + figure; + ax1 = subplot(2, 1, 1); + plot(f12, abs(T12)); + set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude'); + + ax2 = subplot(2, 1, 2); + plot(f12, mod(180+180/pi*phase(T12), 360)-180); + set(gca, 'xscale', 'log'); + ylim([-180, 180]); + yticks([-180, -90, 0, 90, 180]); + xlabel('Frequency [Hz]'); ylabel('Phase'); + + linkaxes([ax1,ax2],'x'); + xlim([2, 500]); +#+end_src + +#+NAME: fig:tf_geophones +#+HEADER: :tangle no :exports results :results value raw replace :noweb yes +#+begin_src matlab :var filepath="figs/tf_geophones.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:tf_geophones +#+CAPTION: Estimated transfer function between the two geophones +#+RESULTS: fig:tf_geophones +[[file:figs/tf_geophones.png]]