Sheila Peacock
Dear All,
I had some odd results from the new subroutine "fap.c" for deconvolving
FAP files, and tracked these down to the absence of extrapolation at
the low- and high-frequency ends of the FAP table.
Two postscripts attached: both generated by:
funcgen impulse npts 10000 delta 0.01
transfer to fap subtype <filename>
The FAP files are at the end. The first was for a long-period (LP)
seismometer, the second for a broadband (BB) seismometer.
What I found was that the lower trace in each postscript, which is ugly
with either a spike or a lot of zig-zagging Gibbs' phenomenon, was
generated by the existing fap.c code because it does not extrapolate
the amplitude and phase values beyond the lowest and highest values
provided in the input table, but just leaves them flat at those
values down to zero frequency and up to the Nyquist. For the first
FAP table (the LP seismo), the highest frequency in the FAP table
is 0.3 Hz so all the amplitudes for frequencies between 0.3 and
50 Hz (=Nyquist) are constant at the value for 0.3 Hz. This flat
spectrum translates into a spike, as seen in the lower trace of
the pair in the first file. For the second FAP table (the BB seismo),
the FAP table ends at 10 Hz when the amplitude (16.2) is only slightly
below the maximum amplitude in the table (19.8), at 6 Hz. So for
frequencies between 10 and 50 Hz (=Nyquist) the fap.c code is
using the value at 10 Hz, which is far from zero and is not tapered
to zero at 50 Hz, hence the zig-zaggy ringing.
The upper trace in each pair were derived from my slightly
altered version of fap.c in which I made the code
extrapolate the amplitude and phase backwards from the first
two values in the FAP table (for frequencies below the lowest
frequency in the FAP table), and forwards from the final two
values in the FAP table (for frequencies between the highest
frequency in the FAP table and the Nyquist).
I know this is risky because it depends on the first
two and last two amplitudes in the FAP table establishing a slope
that sends the extrapolated values down to zero rapidly. Also, it
could be argued that we should not be trying to use "transfer" with
a FAP table for frequencies outside the range of the table. Also I
extrapolated both amplitude and phase, but perhaps I should
extrapolate amplitude only, or extrapolate the complex values
instead.
Regards,
Sheila Peacock
AWE Blacknest.
First one (LP seismo)
1.00000E-02 0.394768 282.522
2.00000E-02 3.11046 146.286
3.00000E-02 6.05993 53.0681
4.00000E-02 7.18071 -21.5560
5.00000E-02 6.35654 -82.4200
6.00000E-02 4.83388 -130.514
7.00000E-02 3.48233 -168.334
8.00000E-02 2.49664 -198.765
9.00000E-02 1.81652 -224.067
1.00000E-01 1.34931 -245.792
0.120000 0.792365 -282.307
0.140000 0.499787 -313.204
0.160000 0.333557 -340.907
0.180000 0.232700 -366.774
0.200000 0.168101 -391.681
0.220000 0.124745 -416.286
0.240000 9.43200E-02 -441.157
0.260000 7.18980E-02 -466.803
0.280000 5.44240E-02 -493.556
0.300000 4.01190E-02 -521.284
Second one (BB seismo)
0.0100 7.13173e-05 -135.464
0.0167 0.000240774 -149.642
0.0200 0.000360009 -154.098
0.0208 0.00039232 -155.002
0.0213 0.00041319 -155.539
0.0217 0.000430258 -155.954
0.0222 0.000452058 -156.457
0.0227 0.000474375 -156.941
0.0233 0.000501834 -157.499
0.0238 0.000525282 -157.947
0.0244 0.000554097 -158.463
0.0250 0.00058365 -158.959
0.0256 0.00061394 -159.436
0.0263 0.000650208 -159.969
0.0270 0.000687477 -160.478
0.0278 0.000731292 -161.033
0.0286 0.00077641 -161.562
0.0294 0.00082283 -162.065
0.0303 0.000876605 -162.605
0.0313 0.000938283 -163.172
0.0323 0.00100199 -163.709
0.0333 0.00106772 -164.218
0.0345 0.00114926 -164.794
0.0357 0.00123372 -165.336
0.0370 0.0013285 -165.888
0.0385 0.0014421 -166.485
0.0400 0.00156024 -167.042
0.0417 0.00169962 -167.631
0.0435 0.00185355 -168.21
0.0455 0.00203223 -168.806
0.0476 0.00222852 -169.384
0.0500 0.00246372 -169.993
0.0526 0.0027316 -170.598
0.0556 0.00305759 -171.236
0.0588 0.00342527 -171.854
0.0625 0.00387607 -172.501
0.0667 0.00442114 -173.163
0.0714 0.00507315 -173.828
0.0769 0.00589251 -174.521
0.0833 0.00692249 -175.237
0.0909 0.00825248 -175.985
0.1000 0.00999766 -176.77
0.1260 0.0159006 -178.567
0.2000 0.0401326 177.966
0.3684 0.13623 172.558
0.5000 0.250858 168.934
0.6840 0.469041 164.116
0.8000 0.64112 161.14
0.8618 0.743638 159.565
0.9283 0.862333 157.872
1.0000 1.00000 156.052
1.1447 1.30823 152.381
1.3104 1.71054 148.172
1.4286 2.0292 145.16
1.5000 2.23426 143.334
1.6510 2.6985 139.46
1.8171 3.25583 135.168
2.0000 3.92373 130.406
2.1544 4.52934 126.349
2.3208 5.22192 121.942
2.5000 6.00992 117.149
2.6566 6.73025 112.922
2.8231 7.52378 108.392
2.8571 7.68882 107.462
3.0000 8.39191 103.539
3.3019 9.91157 95.1709
3.6342 11.5957 85.8732
4.0000 13.3927 75.6033
4.3089 14.8058 66.9726
4.6416 16.1685 57.8003
5.0000 17.4074 48.1589
5.3133 18.2741 40.0068
5.6462 18.969 31.6848
6.0000 19.4637 23.2734
6.3164 19.7134 16.1514
6.6494 19.8051 9.0688
6.6667 19.8058 8.713
7.0000 19.7439 2.0663
7.3186 19.5746 -3.9081
7.6517 19.3067 -9.7789
8.0000 18.9506 -15.5323
8.6177 18.1893 -24.8586
9.2832 17.2663 -33.8022
10.000 16.2303 -42.3532
I had some odd results from the new subroutine "fap.c" for deconvolving
FAP files, and tracked these down to the absence of extrapolation at
the low- and high-frequency ends of the FAP table.
Two postscripts attached: both generated by:
funcgen impulse npts 10000 delta 0.01
transfer to fap subtype <filename>
The FAP files are at the end. The first was for a long-period (LP)
seismometer, the second for a broadband (BB) seismometer.
What I found was that the lower trace in each postscript, which is ugly
with either a spike or a lot of zig-zagging Gibbs' phenomenon, was
generated by the existing fap.c code because it does not extrapolate
the amplitude and phase values beyond the lowest and highest values
provided in the input table, but just leaves them flat at those
values down to zero frequency and up to the Nyquist. For the first
FAP table (the LP seismo), the highest frequency in the FAP table
is 0.3 Hz so all the amplitudes for frequencies between 0.3 and
50 Hz (=Nyquist) are constant at the value for 0.3 Hz. This flat
spectrum translates into a spike, as seen in the lower trace of
the pair in the first file. For the second FAP table (the BB seismo),
the FAP table ends at 10 Hz when the amplitude (16.2) is only slightly
below the maximum amplitude in the table (19.8), at 6 Hz. So for
frequencies between 10 and 50 Hz (=Nyquist) the fap.c code is
using the value at 10 Hz, which is far from zero and is not tapered
to zero at 50 Hz, hence the zig-zaggy ringing.
The upper trace in each pair were derived from my slightly
altered version of fap.c in which I made the code
extrapolate the amplitude and phase backwards from the first
two values in the FAP table (for frequencies below the lowest
frequency in the FAP table), and forwards from the final two
values in the FAP table (for frequencies between the highest
frequency in the FAP table and the Nyquist).
I know this is risky because it depends on the first
two and last two amplitudes in the FAP table establishing a slope
that sends the extrapolated values down to zero rapidly. Also, it
could be argued that we should not be trying to use "transfer" with
a FAP table for frequencies outside the range of the table. Also I
extrapolated both amplitude and phase, but perhaps I should
extrapolate amplitude only, or extrapolate the complex values
instead.
Regards,
Sheila Peacock
AWE Blacknest.
First one (LP seismo)
1.00000E-02 0.394768 282.522
2.00000E-02 3.11046 146.286
3.00000E-02 6.05993 53.0681
4.00000E-02 7.18071 -21.5560
5.00000E-02 6.35654 -82.4200
6.00000E-02 4.83388 -130.514
7.00000E-02 3.48233 -168.334
8.00000E-02 2.49664 -198.765
9.00000E-02 1.81652 -224.067
1.00000E-01 1.34931 -245.792
0.120000 0.792365 -282.307
0.140000 0.499787 -313.204
0.160000 0.333557 -340.907
0.180000 0.232700 -366.774
0.200000 0.168101 -391.681
0.220000 0.124745 -416.286
0.240000 9.43200E-02 -441.157
0.260000 7.18980E-02 -466.803
0.280000 5.44240E-02 -493.556
0.300000 4.01190E-02 -521.284
Second one (BB seismo)
0.0100 7.13173e-05 -135.464
0.0167 0.000240774 -149.642
0.0200 0.000360009 -154.098
0.0208 0.00039232 -155.002
0.0213 0.00041319 -155.539
0.0217 0.000430258 -155.954
0.0222 0.000452058 -156.457
0.0227 0.000474375 -156.941
0.0233 0.000501834 -157.499
0.0238 0.000525282 -157.947
0.0244 0.000554097 -158.463
0.0250 0.00058365 -158.959
0.0256 0.00061394 -159.436
0.0263 0.000650208 -159.969
0.0270 0.000687477 -160.478
0.0278 0.000731292 -161.033
0.0286 0.00077641 -161.562
0.0294 0.00082283 -162.065
0.0303 0.000876605 -162.605
0.0313 0.000938283 -163.172
0.0323 0.00100199 -163.709
0.0333 0.00106772 -164.218
0.0345 0.00114926 -164.794
0.0357 0.00123372 -165.336
0.0370 0.0013285 -165.888
0.0385 0.0014421 -166.485
0.0400 0.00156024 -167.042
0.0417 0.00169962 -167.631
0.0435 0.00185355 -168.21
0.0455 0.00203223 -168.806
0.0476 0.00222852 -169.384
0.0500 0.00246372 -169.993
0.0526 0.0027316 -170.598
0.0556 0.00305759 -171.236
0.0588 0.00342527 -171.854
0.0625 0.00387607 -172.501
0.0667 0.00442114 -173.163
0.0714 0.00507315 -173.828
0.0769 0.00589251 -174.521
0.0833 0.00692249 -175.237
0.0909 0.00825248 -175.985
0.1000 0.00999766 -176.77
0.1260 0.0159006 -178.567
0.2000 0.0401326 177.966
0.3684 0.13623 172.558
0.5000 0.250858 168.934
0.6840 0.469041 164.116
0.8000 0.64112 161.14
0.8618 0.743638 159.565
0.9283 0.862333 157.872
1.0000 1.00000 156.052
1.1447 1.30823 152.381
1.3104 1.71054 148.172
1.4286 2.0292 145.16
1.5000 2.23426 143.334
1.6510 2.6985 139.46
1.8171 3.25583 135.168
2.0000 3.92373 130.406
2.1544 4.52934 126.349
2.3208 5.22192 121.942
2.5000 6.00992 117.149
2.6566 6.73025 112.922
2.8231 7.52378 108.392
2.8571 7.68882 107.462
3.0000 8.39191 103.539
3.3019 9.91157 95.1709
3.6342 11.5957 85.8732
4.0000 13.3927 75.6033
4.3089 14.8058 66.9726
4.6416 16.1685 57.8003
5.0000 17.4074 48.1589
5.3133 18.2741 40.0068
5.6462 18.969 31.6848
6.0000 19.4637 23.2734
6.3164 19.7134 16.1514
6.6494 19.8051 9.0688
6.6667 19.8058 8.713
7.0000 19.7439 2.0663
7.3186 19.5746 -3.9081
7.6517 19.3067 -9.7789
8.0000 18.9506 -15.5323
8.6177 18.1893 -24.8586
9.2832 17.2663 -33.8022
10.000 16.2303 -42.3532
Attachments
-
Sheila,
I guess there are a couple of options here. I think the easiest is
the use the freqlimits option in transfer to limit where the transfer
is applied. The other possibility is to implement an option for the
FAP transfer to do what you have done and extrapolate the FAP outside
it's prescribed range.
I think Arthur had some rationale for making the decision we did about
the FAP code and not extrapolating, but I cannot remember what it was.
Brian
On Jun 23, 2011, at 8:13 AM, Sheila Peacock wrote:
Dear All,
I had some odd results from the new subroutine "fap.c" for
deconvolving
FAP files, and tracked these down to the absence of extrapolation at
the low- and high-frequency ends of the FAP table.
Two postscripts attached: both generated by:
funcgen impulse npts 10000 delta 0.01
transfer to fap subtype <filename>
The FAP files are at the end. The first was for a long-period (LP)
seismometer, the second for a broadband (BB) seismometer.
What I found was that the lower trace in each postscript, which is
ugly
with either a spike or a lot of zig-zagging Gibbs' phenomenon, was
generated by the existing fap.c code because it does not extrapolate
the amplitude and phase values beyond the lowest and highest values
provided in the input table, but just leaves them flat at those
values down to zero frequency and up to the Nyquist. For the first
FAP table (the LP seismo), the highest frequency in the FAP table
is 0.3 Hz so all the amplitudes for frequencies between 0.3 and
50 Hz (=Nyquist) are constant at the value for 0.3 Hz. This flat
spectrum translates into a spike, as seen in the lower trace of
the pair in the first file. For the second FAP table (the BB seismo),
the FAP table ends at 10 Hz when the amplitude (16.2) is only slightly
below the maximum amplitude in the table (19.8), at 6 Hz. So for
frequencies between 10 and 50 Hz (=Nyquist) the fap.c code is
using the value at 10 Hz, which is far from zero and is not tapered
to zero at 50 Hz, hence the zig-zaggy ringing.
The upper trace in each pair were derived from my slightly
altered version of fap.c in which I made the code
extrapolate the amplitude and phase backwards from the first
two values in the FAP table (for frequencies below the lowest
frequency in the FAP table), and forwards from the final two
values in the FAP table (for frequencies between the highest
frequency in the FAP table and the Nyquist).
I know this is risky because it depends on the first
two and last two amplitudes in the FAP table establishing a slope
that sends the extrapolated values down to zero rapidly. Also, it
could be argued that we should not be trying to use "transfer" with
a FAP table for frequencies outside the range of the table. Also I
extrapolated both amplitude and phase, but perhaps I should
extrapolate amplitude only, or extrapolate the complex values
instead.
Regards,
Sheila Peacock
AWE Blacknest.
First one (LP seismo)
1.00000E-02 0.394768 282.522
2.00000E-02 3.11046 146.286
3.00000E-02 6.05993 53.0681
4.00000E-02 7.18071 -21.5560
5.00000E-02 6.35654 -82.4200
6.00000E-02 4.83388 -130.514
7.00000E-02 3.48233 -168.334
8.00000E-02 2.49664 -198.765
9.00000E-02 1.81652 -224.067
1.00000E-01 1.34931 -245.792
0.120000 0.792365 -282.307
0.140000 0.499787 -313.204
0.160000 0.333557 -340.907
0.180000 0.232700 -366.774
0.200000 0.168101 -391.681
0.220000 0.124745 -416.286
0.240000 9.43200E-02 -441.157
0.260000 7.18980E-02 -466.803
0.280000 5.44240E-02 -493.556
0.300000 4.01190E-02 -521.284
Second one (BB seismo)
0.0100 7.13173e-05 -135.464
0.0167 0.000240774 -149.642
0.0200 0.000360009 -154.098
0.0208 0.00039232 -155.002
0.0213 0.00041319 -155.539
0.0217 0.000430258 -155.954
0.0222 0.000452058 -156.457
0.0227 0.000474375 -156.941
0.0233 0.000501834 -157.499
0.0238 0.000525282 -157.947
0.0244 0.000554097 -158.463
0.0250 0.00058365 -158.959
0.0256 0.00061394 -159.436
0.0263 0.000650208 -159.969
0.0270 0.000687477 -160.478
0.0278 0.000731292 -161.033
0.0286 0.00077641 -161.562
0.0294 0.00082283 -162.065
0.0303 0.000876605 -162.605
0.0313 0.000938283 -163.172
0.0323 0.00100199 -163.709
0.0333 0.00106772 -164.218
0.0345 0.00114926 -164.794
0.0357 0.00123372 -165.336
0.0370 0.0013285 -165.888
0.0385 0.0014421 -166.485
0.0400 0.00156024 -167.042
0.0417 0.00169962 -167.631
0.0435 0.00185355 -168.21
0.0455 0.00203223 -168.806
0.0476 0.00222852 -169.384
0.0500 0.00246372 -169.993
0.0526 0.0027316 -170.598
0.0556 0.00305759 -171.236
0.0588 0.00342527 -171.854
0.0625 0.00387607 -172.501
0.0667 0.00442114 -173.163
0.0714 0.00507315 -173.828
0.0769 0.00589251 -174.521
0.0833 0.00692249 -175.237
0.0909 0.00825248 -175.985
0.1000 0.00999766 -176.77
0.1260 0.0159006 -178.567
0.2000 0.0401326 177.966
0.3684 0.13623 172.558
0.5000 0.250858 168.934
0.6840 0.469041 164.116
0.8000 0.64112 161.14
0.8618 0.743638 159.565
0.9283 0.862333 157.872
1.0000 1.00000 156.052
1.1447 1.30823 152.381
1.3104 1.71054 148.172
1.4286 2.0292 145.16
1.5000 2.23426 143.334
1.6510 2.6985 139.46
1.8171 3.25583 135.168
2.0000 3.92373 130.406
2.1544 4.52934 126.349
2.3208 5.22192 121.942
2.5000 6.00992 117.149
2.6566 6.73025 112.922
2.8231 7.52378 108.392
2.8571 7.68882 107.462
3.0000 8.39191 103.539
3.3019 9.91157 95.1709
3.6342 11.5957 85.8732
4.0000 13.3927 75.6033
4.3089 14.8058 66.9726
4.6416 16.1685 57.8003
5.0000 17.4074 48.1589
5.3133 18.2741 40.0068
5.6462 18.969 31.6848
6.0000 19.4637 23.2734
6.3164 19.7134 16.1514
6.6494 19.8051 9.0688
6.6667 19.8058 8.713
7.0000 19.7439 2.0663
7.3186 19.5746 -3.9081
7.6517 19.3067 -9.7789
8.0000 18.9506 -15.5323
8.6177 18.1893 -24.8586
9.2832 17.2663 -33.8022
10.000 16.2303 -42.3532
<
lpcompareFAPwithslope
.ps
<
bbcompareFAPwithslope
.ps>_______________________________________________
sac-dev mailing list
sac-dev<at>iris.washington.edu
http://www.iris.washington.edu/mailman/listinfo/sac-dev
-
Sheila and others,
When writing fap.c, we decided not to do any extrapolation beyond the
given limits. The idea is that those would be handled by other arguments
on the TRANSFER line -- just as when doing a polezero deconvolution. (See
the TRANSFER help file at
http://www.iris.edu/software/sac/commands/transfer.html.)
Applying an LP response directly to an impulse is asking for trouble.
Abrupt changes as one sees in that waveform mean very high frequencies
that must be gotten rid of before applying the response. When I did the
following:
SAC> fg impulse
SAC> fft
DC level after DFT is 1
SAC> plotsp am
the amplitude was constant for all frequencies. One cannot apply an
instrument correction to a signal that does not go to zero at the
Nyquist.
I have also worked on the interpolation routine in SAC, and the methods of
interpolation do not work well with such abrupt changes. SAC routines are
meant to be applied only to real data.
If I have misunderstood your points, please let me know.
Arthur
Brian Savage
Arthur Snoke
