- Email: [email protected]
The Fourier transform of a voltammetric peak and its use in resolution enhancement
371
J. Electroanal. Chem., 2% (1990) 371-394 Elsevier Sequoia S.A., Lausanne
The Fourier transform of a voltammetric peak and its use in resolution enhancement * Sten 0. Engblom Laboratory of Analytical
Chemistry, A*bo Akademi,
Biskopsgatan
8, SF-20500 Abe (Finland)
(Received 18 May 1990)
Abstract
The theory presented shows how the Fourier transform of a voltammetric peak can be utilized for finding peak parameters and enhancing the resolution. It is based on the assumption that the voltammettic or polarographic peak can be described by the se&-function. Potential uses of the Fourier transform of this function are discussed, but the main theme is resolution enhancement. New deconvolution-convolution functions for this purpose are introduced. It is possible to retain either the height of the original peak, or the area under it, by choosing the correct deconvolution function.
INTRODUCTION
One of the reasons, besides improved sensitivity, for introducing polarographic techniques with peak-shaped responses was originally to enhance the resolution when several electroactive species were present [l]. The possibilities for a simultaneous determination of several components are limited in classical dc polarography by the absence of a reliable base-line. In ac polarography, and other methods with similar responses, this difficulty is in many cases reduced. Also, these techniques are better suited to handle a situation where the different half-wave potentials are close to each other. However, in some cases serious overlapping of peaks does occur and a simultaneous determination is difficult, if not impossible. Several examples of this situation exist, among them solutions containing lead and thallium, or cadmium and indium. Many efforts have been made to either prevent the overlapping by chemical methods, or to improve the recorded curve using mathematical means. The subject of this article is the manipulation of the recorded data in the frequency domain, i.e. after Fourier transformation. Especially, the possibility of changing the width of the peaks without affecting their area, height or position, is
l
Presented at the J. Heyrovsky Centennial Congress on Polarography, Prague, 20-25
0022-0728/90,‘$03.50
Q 1990 - Elsevier Sequoia S.A.
August 1990.
372
treated. This operation, usually a deconvolution, is one of the major applications of Fourier transformation in ac polarography, according to Smith [2]. The theory presented below applies to reversible ac polarography, but the same type of response occurs in several electroanalytical techniques, e.g. semidifferential electroanalysis [3] and derivative neopolarography [4]. Grabaric et al. [5] introduced this type of resolution enhancement in electroanalytical chemistry and their deconvolution function was later slightly modified by Engblom and Ivaska [6]. However, this early work lacked the proper theory needed to explain the operation in detail. It is the purpose of this work to give that theoretical background, and at the same time suggest some new deconvolution and also convolution functions of use in resolution enhancement. All calculations in this work were performed using the MATLAB program (trademark of The MathWorks Inc.) on a VAX 8800 computer (trademark of Digital Equipment Corporation). The same program was used to generate the figures shown. THE CONTINUOUS
The general
CASE
form of the function
p(x)=cosh2 [b;lx
under
consideration
here is
(1)
- c)]
The main characteristics of this peak-shaped curve, see Fig. lA, are the height a, the width-determining parameter b, and the position c. A usual measure of width is the width at half height of the peak, Wl,2 = ln(3 + \/8)/b. Further, it is easily shown that the area under the curve equals 2a/b. For later reference, please note that a peak with a large b is narrower and the area under the curve is smaller in comparison with a peak with a smaller value for b. The Fourier transform of this function is given by [7] P(y)
=jm
p(x)
exp(-iny)dx=exp(-icy)?
sitiy/tJ2b)
--m Three main qualities of the peak have thus been separated into three factors: The first factor exp( - icy) is a function of the peak’s original position c, the second is the area under the peak, and the third factor is a function of the width of the peak, this time slightly above the half height, W,.,,,,,,. = r/2b. In Fig. 1B transforms are shown for different b’s. Here one can recognize a fact, which one feels intuitively is correct, namely that a sharp peak (large b) has a “wide” transform. Fast changes in a time-domain curve correspond to high frequencies. On the other hand, a very wide peak has a “sharp” transform. Very few, low frequencies are needed to describe the transform of such a curve. The main reason for performing a Fourier transformation is generally to separate signal from noise. Usually, as in this case, the signal is composed of mainly low frequencies and high-frequency noise is well separated. From eqn. (2) it is clear that not only is the noise separated from the signal, but also that some interesting
313
X
1
I
2-
B 1.5 z
l-
0.5 -
,, ...’
0-10 ----.---8
_--6
-4
-2
0
-._ ---_____ -_____ 2
4
6
8
10
Y Fig. 1. Three voltammetric peaks of different widths (A) and their corresponding Fourier transforms (B). The peak parameters were as follows: o = 1 and c = 5 for all peaks, and 6 = 1 ( -),2(------)and 3 (. . . . . .).
374
characteristics of the peak are disentangled. If a system exhibits only one peak it would be a very simple task for, e.g., an automated instrument to find the area under the peak using P(0) = 2a/b, and the position of the peak using cy = n7~ + y # 0, n = . . . , arctan{ -Im[P(y)]/Re[P(y)]), (2n - l)a/2 d cy 6 (2n + l)n/2, Here Im[ ] and -2, -1, 0, 1, 2 )...) the n = 0 case being the most straightforward. Re[ ] denote the function’s imaginary and real parts, respectively. These parameters are usually the only ones required in a routine analysis. Note that a conventional measurement of peak height and position in many cases involves the use of only one point on the curve. A more elaborate method would use some kind of curve fitting in the region of interest. In the procedure described above no fitting is needed and all recorded data points are used. Furthermore, the data used would be essentially noise-free. This approach has been used previously in the analysis of voltammetric waves [S]. However, the main theme of this article is the manipulation of the third factor of the right-hand side in eqn. (2). This factor determines the width of the time-domain peak and is readily accessible in the frequency domain. Multiplication with a suitable function would clearly be the most straightforward way to accomplish a desired change in width. It would thus be possible to affect the width without any knowledge of either a or c. In earlier attempts to change the width of the peak in this way, b was also assumed unknown. This, however, produced some unwanted side-effects. These early methods will be reviewed to expose the problem. No knowledge of any peak parameters available The multiplication of two functions in the frequency domain is equivalent to a convolution operation in the time domain. Thus, if U(y), V(y) and W(y) are the Fourier transforms of U(X), v(x) and w(x), then the convolution U(X) = /@‘,v( z)w(x - z)dz and the inverse transform of U(y) = V( y)W( y) yield the same result. Finding v(x) when U(X) and w(x) are known is called a deconvolution. By necessity (simply because the transform, eqn. (2), was not known), the earlier methods were deconvolutions. Deconvolution functions were defined in the time domain and their transforms were calculated numerically. The transform of the voltammetric peak was divided by the transform of the deconvolution function and, if the form of the deconvolution function had been chosen suitably, a narrower peak was the result of an inverse transform. However, the height of the new peak was not predictable, it was a function of the width of the deconvolution function but the exact relation was not known. The deconvolution functions used were of the same functional form as the peak but the peak height was either arbitrary [5] or equal to N, the number of points recorded [6]. Here, as an example, let us consider a deconvolution function with a peak height equal to b’/2, where b’ is the width-determining parameter in the deconvolution function. We have thus dpa(x)
=
b’/2 cash* ( b’x )
the transform
of which
is given
by eqn. (2) by letting
c = 0 and
a = b’/2,
and
375
replacing b with b’. The name of the function is derived from deconvolution (‘d ‘) and voltammetric peak-shape (‘p’), while the superscript ‘a’ indicates that this particular function does not affect the area under the peak it is operating on. The reason for this is that the area under the curve defined in eqn. (3) is equal to one. Also, this deconvolution function will not affect the position of the peak, because it is centered over the origin. The transform of d@(x) consists, as P(y), of three factors, but now the first two are both equal to one. Capital letters are used to indicate Fourier transforms, so that, e.g., the transform of &P(x) will be denoted gPa( y). The deconvolution is accomplished by, in the frequency domain, forming the quotient 2a b’ sinh(ry/2b’)
RP”(y) =p(Y)/Dp”(Y) =exd-icy)~
b
sinh(1TY/2b)
Only the third factor has been affected by the operation, the first and second factors, determining the position of and area under the peak, are still intact. The inverse transform of the result RP(y) is given by [7]
rpaCx) = ‘iii’
cosh[2b(
sin( oh/b’) x - c)] + cos( ah/b’)
assuming b’ > b. ‘r ’ in the function name refers to “result”. M’=(y) lacks an inverse transform in the case 6’ < b (see the Appendix). The result of a deconvolution using a deconvolution function of the same shape as the voltammetric peak is thus also a peak. In this case the area under it is still the same as under the original peak, 2a/b, but the peak height is now a complicated function of b’, given by a tan(rb/2b’)/(ab/2b’). The width of the peak at half height is also, as could be it to equal (2/b)ln[cos(ab/2b’) expected, a function of b’; one finds + cos’( rb/2b’) + 11. The useful range of b’ has already been limited to b’ > b. It is now of interest to investigate the behaviour of rp”(x) for different values of b’ and especially the limits b’ 4 b and b’ + 00. As b’ + b it is evident that the width of the peak tends to zero while the height tends to infinity, while the area remains constant. Thus, we have rp”(x) -+ (2a/b)6(x - c ) as b’ + b, where 6(x) is Dirac’s delta function. The other limit is likewise interesting as one can show that rp”(x) --, p(x) as b’ + 00. The resulting peak will thus vary between a delta function, which must be seen as the ideal response for a resolution-enhancing operation, and the original peak itself, as b’ varies in the interval [b, cc). Of course, once the analytical form of P(y) is known, the possibility of using convolution instead of deconvolution, is opened up. In the deconvolution operation the transform of the deconvolution function appears in the denominator in the quotient shown in eqn. (4). Since DPa(y) is itself peak-shaped and decreases rapidly towards zero as 1y 1 increases, there can be difficulties in implementing this operation in practice. A convolution, i.e. a multiplication, using c-P”(y)
=
sinh( ry/2b’) vy/2b’
376
X
The two component peaks are also shown Fig. 2. A curve with two similar, overlapping peaks ( -). (- - - - - -). The peak parameters were: a = 1 and b = 1 for both peaks, and c = 4 for the first peak and c = 6 for the second.
would on the other hand decrease the risk for the noise-enhancement inherent in the division. However, since CP”(y) + cc as y + cc it is clear that this operation will also enhance noise if carried out without proper precautions. More will be said on these problems in a later section. There is no time-domain equivalent to CP(y) (see the Appendix). A system consisting of several peaks, see e.g. Fig. 2, does not present a problem if one assumes that the different responses are additive. That is, if a curve consisting of several peaks can be written as
mp(x) =Plb> +p*(x) + ... +P”(x> where each pi(x), i = 1,2,. . ., n, is given by eqn. (1) after inserting parameters cli, bi and ci. The linear property of the Fourier transformation ensures that the transform of mp( x) is given by
where each P,(y) is given by eqn. (2) by inserting the proper parameters. Now, a multiplication with U”(y) or a division with DP(y) will, simply because of the distributive nature of those operations, affect each peak individually. In the inverse transform the effect on the peaks is determined by their original width, in accordance with the relations given above (for an example of a situation with two similar
311
1.8,
X
Fig. 3. Theoretical results of a resolution enhancement of the curve shown in Fig. 2. The use of eqn. (3) will preserve the area under the peak (A) while eqn. (9) leaves the height unchanged (B). The curves shown are: The original curve (), the result using b’ = 2 (- - - - - -) and the result using b’ = 1.5 ( . . . . ..).
378
peaks see Fig. 3A). If complete separation which is unchanged, is a useful measure.
is accomplished
the area under
the peak,
The width-determining parameter b assumed known In reversible voltammetry knowing b corresponds to knowing n, the number of electrons participating in the electrode reaction. Note that in a one-peak system b is relatively easy to find using the third factor in eqn. (2). One could, e.g., use the relation x = In (x/y + 7 1 + (x/y) ), which is valid for y = x/sinh(x), and find x in an iterative manner. This way of determining b would benefit from the advantages mentioned earlier, all points on a curve would contribute and the influence of noise would be eliminated to a large extent. Once b is known, it is a simple matter to introduce a factor into the convolution-deconvolution functions to achieve constant peak height in the result. In this case, dp becomes dph(x)
=
Similarly,
CP”( y)
tan( nb/2b’) ?rb/2b’
b’/2
the convolution =
mb/2b’ tan( nb/2b’)
(9)
cash* (b’x) function
is given by
sinh( vy/2b’) vy,‘2 b’
The superscript ‘h’ indicates b’. The resulting peak
(10)
here that the peak height of the result is independent
sin( nb/b’) rph(x)
= tan(ri/2b’)
cosh[2b(x
- c)] + cos(rb/b’)
of
(11)
is, of course, of the same shape as rp”(x) but the peak height is now equal to a for all b’ while the area under the peak is a strong function of b’; an integration yields the result ma/b’ tan(ab/2b’). The width of the peak at half height is still the same as that of rp”(x), or (2/b) In [cos(sb/2b’) + cos*(ab/2b’) + 11. rph(x) behaves in the same way as rp”(x) as b’ becomes very large, i.e. we have that rph(x) + p(x) as b’ + co. The situation is different at the other limit. As b’ -B b we see that both the width of the peak as well as the area under it tend to zero, but the peak height remains constant. This could lead one to expect a limit consisting of a finite delta function at x = c. That this limit in fact cannot be reached is evident in the frequency domain where we have RPh(y)
= exp(-icy)
b’ sinh( Ty/2 b' ) b’ tan(yb/2b’>
b sinh( ry/2b)
(12)
and it is clear that RPh(y) + 0 for all y as b’ + b. A multi-peak system is more problematic in this case. In Fig. 3B the simple case of two similar peaks is used, but, unfortunately, it is highly unlikely that all peaks have the same width. One then has to choose one b and accept that only peaks of
319
this width will retain their true heights. As an example, consider a system consisting of n peaks. Each peak p,(x) is characterized by the parameters ai, b, and ci, i=12 7 ,.*., n. If b, is chosen to be used in the resolution-enhancing operation, we get up(x)
=
tan( nbj/2b’)
b’/2
nb,/2b’
cosh’( b'x )
(13)
or CP”( y) =
rbj/2 b’ tan( sbj/2b’)
If the effect on peak p;(x) r#*(X)
=
sinh( vry/26’) I ry/2b i #j, one arrives at the result
is studied,
abj b, tan( abj/2b’)
(14)
sin( n-b/b’) (15)
- c)] + cos(rb;/b’)
cosh[2b,(x
b, tan( mbJ2b’) = bi tan(rbj/2b’) The factor relating
(16)
Vh(4
rph*(x) and rph(x) /> 1
if bi > b,
= 1
if bi = b,
bj tan( ab/2b’)
< 1
if bi < b,
bi tm( rbj/2b’)
’ -B I --+O +oo
in eqn. (16) has the following
as b’ + 00 as b’-+b,, b,
properties
(17)
bi>b,
In the last two cases it is assumed that it is meaningful to let b’ approach either bi or b,. In a system consisting of several peaks b’ should always be greater than, or equal to, the largest of the bi’s. This is required in order to ensure that inverse transforms of all the convolved peaks really exist (see the Appendix). Thus, in a multi-peak system, deconvolution using dph(x) or convolution using CPh( y) will produce the correct result only if the b in the deconvolution-convolution function is the same as in the original peak. The resulting peak will be smaller if the b used is larger than the b of the original peak, and larger if b in the deconvolution-convolution function is smaller than that of the peak. A natural question at this stage is if there are any other possible deconvolutionconvolution functions. Consider eqn. (2) again. It is evident that if b is known it is a simple matter to replace it with b’ while preserving the original shape of the peak. That is, as before it is possible to change the width of the peak without affecting, or needing to know, the area under it (alternatively, the height) or its position. To this one can now add that it is also possible to preserve the functional form of the peak.
380
All one has to do is to multiply b sinh( ry/2
the transform
of the peak, eqn. (2) with the function
b)
08)
cR”( ‘) = b’ sinh( ?ry/2b’) or, equivalently, divide the transform with DRa( y) = l/CP( y). function does not have an inverse transform (see the Appendix). transform of DRa( y) is given by [7] sin( Irb/b’) b’ dracx) = G cosh(2bx) + cos( mb/b’)
The CZV’( y) The inverse
(19)
What we have here is nothing but the result of a deconvolution using u”“(x) shown in eqn. (5), but for a constant factor. Hence the naming of this function, the first convolution letter in the function name, ‘d ‘, ‘c’ or ‘r’, still refer to deconvolution, and result, respectively, while the second letter, ‘r ‘, indicates that this function actually is of the same form as the result of a deconvolution-convolution of “the first kind”. The result of the operation is here
RR”(y)
= P( y)CP(
This function
y) =
has an inverse
P( y)/DP( similar
y) = exp(
--icy)?sinl(/tJib,) (20)
to eqn. (l),
(21) This new peak is of the same shape as the original one; they have the same positions but the new peak is narrower by the factor b/b’. Furthermore, the height of the new peak is increased by the factor b’/b. The area under the peak, however, remains the same (hence the superscript ‘a’). Clearly, it is also possible to derive functions of “the second kind” that will preserve the height of the original peak. These are simply defined by CRh( y) = (b/b’)CF( y) and DRh( y) = l/CRh( y). If P(y) in eqn. (2) is now operated on by either of these two, the result is
RRh(y)
= exp( -icy)%
sitiy/t!..b,)
with the time-domain
equivalent
rrh(X) = cosh2[b;x
- c>j
(22)
(23)
which gives a peak of the same height as the original one, but narrower by the factor b/b’. The area under the peak is also decreased by the same factor. In an operation of “the second kind” b’ is still restricted to the interval [b, 00). A b’ smaller than b would not be forbidden but it would produce a peak wider than the original. Also, in this case the convolution function would have an inverse
381
transform while the deconvolution function would lack one. However, the normal case is a b’ which is larger than b; in fact the larger b’ the narrower will the resulting peak be. As b’ + co we see that rP(x) + (2a/b)6(x - c). The fact that RRh(y) tends to zero for all y in this limit shows that the finite delta function, which is taking the value a for x = c but is zero everywhere else, is not attainable in the time domain. Of course, in both cases the resulting peaks tend to p(x) as
b’ -+ b. In a system exhibiting several peaks, the case where b in the deconvolution-convolution function is incorrect has to be considered once again. In an operation of the second kind this applies to both types of functions, i.e. both the area- and the height-preserving function. Assume that the width-determining parameter of peak j is used while the effect on peak i is studied, i #j. If, as an example, CRa( y) is used, one obtains
2a
RR”*(y) =exp(-icy)7
rY/2bi sinh( ry/2bi)
bj sinh( rry/2 bj ) 6’ sinh( ry/2b’)
bj sinh( Ty/2 bj ) = bi sinh( lry/2b,)
(24)
(25)
RR”( ‘)
The case is even more complicated here than in eqn. (16) since the proportionality factor is a function of y. The main properties are ‘>l
bj sinh( vy/2bj) bi sWry/2bi)
if 6, > b,
= 1
if bi = bj
< 1
if bi < bj
’ + 1
(26)
asi yl -+O asly
+m,
b,
--,CCI aslyl
*co,
bi>bj
-+O
If bi -c b, there is no problem regarding the existence of an inverse transform; we have I RR”*(y) I G I RRa( y) I for all y. This is sufficient since absolute integrability ensures the existence of the inverse (see the Appendix) and /‘?, I RRa*( y) Id y G /?, 1RF(y) ldy = 4nab’/b [9]. However, so far no explicit expression for the time-domain equivalent of RR”*(y) has been found. Something can still be said though about V”*(X). It is the convolution of rP(x) (eqn. 21) and a function given by eqn. (19) if b and 6’ in that equation are replaced by b, and b,, respectively. These two functions are symmetrical and peak-shaped and their convolution yields a symmetrical peak. This peak is centered at x = c as the original peak and the area under it is still 2a/b,. The width of this peak will be determined by not only b’ but also 6,. If b, is very close to bi so that the proportionality factor in eqn. (25) is only slowly decreasing towards zero, it will be possible to get a relatively narrow peak by using a large b’. On the other hand, if the difference in 6, and b, is large, even a very large b’ will produce only a minor effect on the width of the resulting peak.
382
The proportionality factor is thus in a way working like a digital filter where one can affect the cut-off frequency by varying bj. Similar conclusions are valid when considering the use of a height-preserving function. In this case the relation between RRh*(y) and RRh(y) is given by b,’ sinh( ry/2 RRh*( y) =
b, )
b: sinh( ny/2 b, )
The properties of the proportionality eqn. (26), but still different enough of them. We have b,? sinh( ry/2bj)
(27)
RRh(y)
factor in eqn. (27) are similar to those given in to warrant a separate exposition of at least some
> b,/b,
if b, > bj
< b,/b,
if 6, < bj
(28)
b,? sinh( ry/2b,) i -b,/b,
aslyl
-0
The missing ones (the b, = bj case and the limits as ) y 1 -+ co) are the same as in eqn. (26). In the b, < b, case the existence of the inverse can, as above, be drawn from the fact that /‘E’, 1RRh*( y) Id y d (b/b,)/‘?, I RRh( y) Id y = 4nab,/b, [9]. The shape of rrh’(x) is similar to that of rr”/*(x), but the two functions are scaled differently; we have rr”*(x) = (bb,/b’bi)rra*(x). If b, > 6, we have to look more closely at the behaviour of the frequency-domain function. As has been stated several times, absolute integrability guarantees the existence of an inverse. If the function under investigation is even (as RR” *( y ) and RRh*( y) are) the consideration can be restricted to the interval [0, cc). Furthermore, the integral can be divided into two parts, lo” ] g(x) I dx = 10”] g(x) I dx + /p ] g(x) I dx. If g(x) in this example is continuous on the closed interval [O,k] then ] g(x) ] is integrable. It is then the second integral which is of interest. When either RRa*( y) or RRh*( y) is considered, assume a k which is large enough to justify the approximation sinh( x) = exp(x). The integrand is then simplified and its behaviour for large y is given by (the constant factors are dropped and c = 0 is assumed for simplicity)
----7
= ”
i
0
if 6, -c bj
0
if b’ < b,b,/( bj - 6,)
and
bi > bj
00
if b’>b,bj/(bi-b,)
and
b,>b, (29)
It is clear that an integral will converge in the first two cases but not in the last case. If b, > b, the magnitude of 6’ is severly limited, which is a considerable disadvantage when considering that in a resolution-enhancing operation of the second kind a large 6’ means better resolution. In practise, dealing with experimental data one is of course obliged to use a discrete version of the results presented above. This will be the theme of the next section.
383
THE DISCRETE
CASE
Assume that an experimental are given by Pk = cosh2[b(xi:+
kh - c)]
curve is known
k=0,1,2
in N equidistant
points.
N-l
)...)
The points
(30)
where xin is the position of the first measured point and h is the distance between consecutive points. These points are of course related to the continuous function of the last section through pk =p(xi,, + kh). It is assumed here and below that the number of points, N, is large enough to characterize the peak and its Fourier spectrum completely. A specific number is discussed below. The discrete Fourier transform is defined by N-l
PI = c
pk exp( - i2nkl/N)
i=O,
1, 2 ,...,
N-
1
(31)
k=O
and the original
data are recovered
through
the inverse
transform
N-l
1 pk = x c P, exp(i2rkl/N) I=0
k=O,
1,2 ,...,
N-l
(32)
The discrete Fourier transform of the array of data in eqn. (30) should related to P(y) in eqn. (2). A reasonable model is P, = exp( - i2rQ/N)
7~‘l/bhN sinh( r’l/bhN
)
I=O,
1,2 ,...,
[N/2]
be closely
(33)
In the discrete case the position of the peak is taken relative to the first sample, i.e. if we have c = xi,, + .$h, 0 G [ G N - 1, then the first factor is a function of E instead of c. The second factor is now the sum of the recorded points rather than the area under the curve, while the third factor is the same as in eqn. (2) if y is replaced by 2rl/hN. It is known that, since pk is a real array, PN_, = P,* (the complex conjugate of P,) [lo]. The points given in eqn. (33) are then enough to generate the whole array. In eqn. (33) [N/2] signifies the greatest integer less than or equal to N/2, which is N/2 if N is even and (N - 1)/2 if N is odd. When a numerically calculated transform is compared with eqn. (33) the fit is perfect except for the highest frequencies, see Fig. 4. The reason for this deviation is not yet known. Equation (33) can now be used to estimate the minimum number of samples required. The third factor in eqn. (33) is rapidly decreasing towards zero. Only P, of small values on I are needed to describe the peak in the Fourier domain. Assume that PI with I> A does not contribute significantly to the transform, and that 1P,, 1 ( O.OOOlP,. The problem is then to find those values of z that satisfy z/sinh(z) < 0.0001. The answer is z >, 12.4231. In eqn. (33) z is r21/bhN. The product hN can be seen as the interval in x used when sampling the function and let us assume this is equal to five peak-widths, or hN = 5 ln(3 + 6)/b. This gives
384
-21 0
20
40
60
80
100
120
20
40
60
80
100
120
-l6ot.
v Fig. plot vs. h=
and eqn. (33) (....e.). In a normal 4. Comparison of a numerically calculated transform ( -) of 1P, 1 vs. I the fit seems perfect (A), while a small discrepancy shows up if, e.g., ln( 1P, I) is plotted I (B). The peak parameters were: a = 1, b = 3 and c = 5. The peak was sampled using xin = 0, 0.078125 and N = 128. The peak is shown in Fig. 1 (. . . . . .).
385
I = 62.1155 ln(3 + a)/;* = 11.0941. The same number of data is needed for the above-N/2 part of the array, so this crude analysis gives a lower limit for the number of samples at about N = 25. Of course, in order to sample high-frequency noise correctly in accordance with the sampling theorem, a much higher sampling rate is usually required. As in the continuous case, it is now possible to find E, and thus c, by using the relation 2a&/N = n?~ + arctan { -Im[P,]/Re[ P,]} which is valid when (2n 1) N/4 Q 51~ (2n + 1) N/4, I f 0, n = . . . , - 2, - 1, 0, 1, 2,. . . . Once again, n = 0 simplifies the calculation even though it restricts the magnitude of (1 considerably. Actually, even if only PI is used, 5 is still restricted to the first quarter of N. This restriction can be eased by adding zeros to the end of the data array before Fourier transformation. PO is useful when an estimate of the area under the peak is required. We have that P,, = If:: pk and therefore hP,, = 2a/b. Here it is necessary that the peak is entirely within the sampled region (meaning that the peak’s contribution to samples outside this region is negligible). It is also possible to find b as explained earlier. Resolution enhancement is achieved in the same manner and using essentially the same functions as in the previous section. A convolution is, in the discrete case, defined as a sum, uk = Cf:i u, w,_,, and the equivalent operation in the frequency domain is a term-by-term multiplication of the discrete Fourier transforms, U, = YW, [ll]. A deconvolution is then defined as a term-by-term division of two Fourier transforms. If a time-domain deconvolution function is used, one has to consider the periodicity implicit in eqns. (31) and (32). Both pk and P, are seen as periodic and, in order to avoid this periodic repetition being “seen” in a arrays, convolution-deconvolution operation, one has to destroy the periodicity. This is accomplished by zero-filling. An aperiodic array is created by appending N zeros to the original array [ll]. In eqn. (3) a deconvolution function was given that did not affect the area under the peak nor was any knowledge of the peak parameters required. In the discrete case this function has to be slightly modified in order to produce the same end result. The discrete transform is, as the continuous one, a product of three factors, see eqn. (33). In the transform of the deconvolution function it is necessary for the first two factors to equal one. The first factor has to equal one, i.e. 5 = 0, in order not to affect the original peak position. The maximum of the deconvolution function has then to occur at x = xi,, instead of at x = 0 as in the continuous case. The second factor, the sum of all points, is related to the area under the curve. If this sum is multiplied by h one gets the Riemann sum, which, if h is small enough, gives a good approximation to the area under the curve. For the deconvolution function not to affect the sum, and thereby the area under the peak, its own sum has to equal one. We can now define dp,” =
hb’/2 K=[-;+I]
,...,
-l,o,l,...,
I;]
cosh*( b’Kh)
and by using
dpY, = dpi_,
(this can be deduced
from eqn. 32) one can construct
386
must be done to the the complete array d’i, k = 0, 1, 2,. . . , N - 1. Zero-filling array in eqn. (34) before rearrangement. The relation to the continuous function is dp,” = h dp”(x, + oh). The area under the curve h dp”(x) is equal to h, and therefore, approximating the area under dpz with hC, dpi, it follows that C, dp,” = 1. The Fourier transform of dpl is then given by DP,a =
r21/b’hN sinh( r’l/b’hN)
1=0,1,2
,...,
[N/2]
(35)
with the reservation, as mentioned earlier, that the fit at the highest frequencies is not perfect. The result of a deconvolution using dpt yields a result similar to that in the continuous case. We get b’ sinh( a’l/b’hN)
RPP = exp( - i2r@/N)
b sinh( a21/bhN )
l=O,1,2
,...,
[N/2] (36)
or after an inverse 2ab’ rpi = rb cosh[26(
transformation sin( rb/b’) xi” + kh - c)] + cos( ah/b’)
k=O,
1, 2 ,...,
N-l
As in the continuous case, b’ can now vary in the interval [b, co). Equation describes rpt accurately for intermediate values of b’. However, as b’ tends there are deviations from this equation. In the limit we have
(37) (37) to b
if k = .$ (38) otherwise where 6,, is Kronecker’s delta function. In the other limit, b’ + 00, we have as before that rpt --, pk. In Fig. 5 the height of and area under rpi are compared with the height and area of the result in a numerical calculation. The fit is quite good, only a minor difference in peak height for 6’ close to b being detected. The peak height in this region is actually growing a little faster than predicted, before attaining a finite value at b’ = b. Note that a large b’ gives a very sharp dp,“, thus necessitating a large number of samples (small h) in order to achieve a correct transform. This would not be necessary if the operation started off directly with DPF or CP: = l/DPp. Using the present theory this is of course possible only if the highest frequencies do not contribute to the signal in any significant way. The other relevant properties of the peak in eqn. (37) were given in the previous section. It is now possible to examine a little more closely the effect of an early deconvolution function, N/cosh’(b’Kh) [6]. The peak height of this function should, however, be discussed first. In the context where this function was defined, slightly different versions of eqns. (31) and (32) were used. The discrete Fourier transform was defined as P, = (l/N& pk exp( -i2nkl/N) while the inverse was calculated according to pk = C, P, exp(i2nkl/N). The difference lies in the position-
ing of the factor l/N. This might be thought of as an insignificant distinction, but it is of some importance. Let once again U,, V, and W, be the discrete transforms of uk, vk and wk. The equivalent operation to U, = V,W, is then uk = (l/N)C, v,w,_,. If, however, eqns. (31) and (32) are used, it is uk = C, v,w,_,. The presence or absence of the factor l/N determines the peak height of the deconvolution function. Therefore, using the conventions of this article, this early deconvolution function is given by dp; =
1
k=[-;+l]
cosh’( b’kh)
The transform
)...,
)...,
is given by (with the same reservation m21/b’hN sinh( 1~‘l/b’hN
Forming
-1,&I
P,/DP,’
the quotient
RP,’ = exp( - i2n[f/N
f=O,
)
(39)
as earlier) IN/21
(40)
gives
) “;
) (crp;
1,2 ,...,
I;]
sffL[ ~~~~~~)’
I = 0, 1, 2,. . . ,[ N/2]
(41)
k
where Crp: = Cp,/Edpi. rp; =
given by
sin( sb/b’)
a( b’)2h mb
The inverse
cosh[2b(
xi” + kh - c)] + cos( ah/b’)
k=O,
1,2 ,...,
N-l (42)
again applies only for intermediate values of b’. As b’ -+ b we have that rpz + as,*, instead of the limiting infinite delta function indicated by eqn. (42). At the other limit (b’ + 00) we have, for a fixed N (and h), that dp: + Sk,, the transform of which equals unity for all 1. A deconvolution has no effect and the result equals the original. In the case where b is known to the experimenter, it is also, in the discrete case, possible to use deconvolution-convolution functions that will preserve the peak height and/or the functional form of the original peak. The discrete version of eqn. (9) is given by dp, =
tan( rb/2b’) nb/2b’
dp” K
.=I-;+1]
,...,
As before, the range of b’ where this function the use of this function is rpkh= tan(rij26’) k=O,
1,2 ,...,
cosh[2b(xi, N-l
-l,o,l,...
:I;]
can be used is limited.
The result of
sin( nb/b’) + kh - c)] + cos(rb/b’) (44)
18 -
A
16 14 12 . 10 86-’ 42-
* . ‘* l**.
2
l *.*.
3
l
4
. 5
.
*
.
l
.
6
1
8
9
10
“--‘-1.2 -
B
I*- . . . . . .t.
. . . . .
l
.
.
*
l
. i
0.8 -
0.6 -
0.4 -
0.2 -
2
3
4
5
6
7
8
9
10
b’
Fig. 5. Resulting peak height (A) and area under the peak (B) as a function of the b’ used in a deconvolution. The predicted height and area of the result rpi (eqn. 37) (. . . . . .) is compared with the result in an actual numerical calculation ( l). In the numerical example eqn. (30) was used to generate synthetic data. The peak parameters were: (I = 1, 6 = 2, c = 5, xin = 0, h = 0.078125 and N = 128 (256 after zero-filling). Equation (34) was used to construct the deconvolution function. The simple formula hE, rp; wasusedto calculate the areas after deconvolution.
389
1.4
.
A
1.24 l
1-
s
0.8 -
: &
0.6 -
..**.*. .
l
l
. .
.
.
l
.
l
.
l
8
9
0.4 -
0.2 -
01 2
3
4
5
6
7
b’
1.4. B
1.2 c
l. 0.8 -
.
*
.
*
.
l
I
8
9
l
I 0.61
l
;
I .
l
l
0.4 -
.
’
t 0.2 - ’ l
01 2
3
4
5
6
10
b’
Fig. 6. A comparison of peak height (A) and area under the peak (B) with the predicted values (. . . . . .). The height-preserving deconvolution function dp,h in eqn. (43) was used in a numerical calculation and the results ( l) were compared with e-qn. (44). The peak parameters were as in Fig. 5. The areas were calculated using hE, rp*“.
390
1.8 1.6
A
1.4
0.8
I-
O-
I
2
3
4
5
6
I
8
9
10
6
7
8
9
10
were generated according Fig. 7. Results of a deconvolution using drxa(eqn. 45). Synthetic data ( -) to eqn. (30) using D = 1, h = 3, c = 5, h = 0.078125 and N = 128 (256 after zero-filling). Results using b’ = 4 (- - - - - -) and b’ = 5 (. . . . .) are shown (A). A logarithmic analysis (B) shows the correct slopes long linear parts are due to the fact that (3, 4 and 5) and the correct intercept (x = 5). The exceptionally the original peak was computer-generated and noise-free.
391
See Fig. 6 for an example of how well the height of and area under the resulting peak in a real calculation adheres to the height and area of the peak in eqn. (44). The fit is again very good except for a deviation in height when b’ is very close to b. The height is, as before, growing faster until the tangent-factor in the denominator forces the peak height towards zero. Considering a multi-peak system, the discussion in the provision section still applies wherever eqn. (44) is valid. Discrete versions of the functions that preserve the functional form of the peak, are also available. They will be given as deconvolution functions below, but the convolution functions can also be used if one takes into account the deviating behaviour at the highest frequencies (e.g. by using filters). The area-preserving deconvolution function of the second kind is thus given by I drKa= c
sin( vrb/b’) cosh[ 2b( xin + id)]
+ cos( ah/b’)
.=[-;+I] ,...,-1,0,1,...,
[;J
(45)
1.8 peak 1
2
1.61.41.2 I-
“2 c
0.8 0.6 0.4 0.2 0
\t.. 0
1
2
3
4
5
6
I
8
9
I 10
X
Fig. 8. Deconvolution of a double-peak system. The peak parameters were for peak 1: a = 1, b = 2 and c = 4.5, and for peak 2: CI= 1, b = 3, and c = 5.5. The peaks were generated separately using eqn. (30) with h = 0.078125 and N = 128 (256 after zero-filling), and then added. In the deconvolution function b of peak 2 was used, i.e. b = 3. The results of using b’ = 4 (- - - - - -) and b’ = 5 (. . . . . .) are shown together with the original curve (). The resulting peak heights of peak 1 cannot yet be predicted. The peak heights of peak 2 are (the true values are shown within brackets): 1.071 (l.C!C@). 1.369 (1.333) and 1.674 (1.667).
392
and the result by ab’/b “r’=
cosh2[b’(xi,+kh-c)]
k=O,
1, 2 ,...,
N-l
The height-preserving function is simply related to dr,” through drf’ = (b’/b)dr:, and the result is similarly given by ori = (b/b’)rr,“. As before, eqn. (46) is valid for a wide range of b’. However, for very large b’ and for b’ close to b there are deviations. As 6’ + co we have that dr,” --+ (hb/2)/cosh2(b~h). We have thus for large 6’ that rrt is following closely the behaviour of rpi for b’ close to b. rrf tends to zero in this limit as a consenquence of the factor b/b’ and the fact that the limiting delta function is finite. In the less interesting limit, b’ + b, both rr: and rrt tend to pk according to eqn. (46). However, the use of eqn. (45) are not practical in this limit. In Fig. 7 some results of the use of eqn. (45) is shown. A logarithmic analysis [12] of the peaks was performed to show that the curves were of the right shape. Finally, Fig. 8 shows a double-peak system with different b’s and the response to a deconvolution using drKa. b of peak 2 was used in the deconvolution function and it is thus that peak which shows a predicted response to the deconvolution. The resulting peak heights of peak 1 are not yet predictable; they are however still directly proportional to a of that peak. The areas under both peaks remain the same and if complete separation is achieved, these can be used. CONCLUSIONS
In this work the theoretical background for resolution enhancement of voltammetric and polarographic peaks by deconvolution has been given. Utilizing the Fourier transform of a peak, several functions were derived, each with a distinct property distinguishing it from the others. The user can choose between two forms of the resulting peaks, and whether the original peak height or the area under the peak shall be retained. Numerical calculations using synthetic data yielded results that agreed well with the theoretical predictions. The present theory works well both in the continuous and the discrete case. Some refinement is still needed in the discrete case for the high-frequency region, but in practice, when filters are in any case needed to control instrumental and other noise, the theory given here is fully adequate. If compared with other methods of finding the components of a fused system of peaks, the present method has the advantage of not needing any initial guesses with regard to number of peaks, positions of peaks, peak heights, peak widths, etc. No iterations of any sort are needed in the resolution-enhancing operation, thus eliminating cases of lengthy, diverging calculations. Among the disadvantages, one can mention the fact that the present theory concerns only reversible peaks; the incorporation of quasi-reversibility into the theory can be very difficult. In the crucial deconvolution-convolution step a sensitivity to noise is characteristic.
393
APPENDIX
The inverse
transform
of a function
P(y)
is given by (Al)
Generally, the integral in eqn. (Al) exists for all x if P(y) is absolutely integrable, or at least for all x # 0, if P(y) converges monotonically to zero as y + f cc [13]. The functional form of interest here occurs in several of the definitive equations in this article. CRa( y) (eqn. 18) is chosen as an example, but the results below apply to all cases after obvious modifications. That this function is absolutely integrable when b’ < b is well known; the integral can be found in most handbooks, see e.g. ref. 14,
If b’ > b, the integrand above will grow beyond any limits as y + cc. To see this use, e.g., the approximation sinh(x) = exp(x) which applies for very large x. Thus, the function is absolutely integrable and an inverse will exist if and only if b’ < b. The inverse is in fact given by eqn. (19) by letting b and b’ change places. The reasoning above could equally well be applied to RPa( y ) or DRa( y ) and explain why the time-domain functions rp”(x) and d?(x) exist when b’ > b but do not when b’ -c b. The functions CRh( y), DRh( y) and RPh( y) differ from the above-mentioned functions only by a constant factor, and are therefore also subject to the same conclusions regarding the existence of their time-domain equivalents.
REFERENCES
1 G.C. Barker in P. Zuman (Ed.), Progress in Polarography. Vol. 2, Interscience, New York, 1962, Ch. 19. 2 D.E. Smith in A.G. Marshall (Ed.), Fourier, Hadamard, and Hilbert Transforms in Chemistry, Plenum Press, New York, 1982, Ch. 15. M. Goto and D. Ishii, J. Electroanal. Chem., 61 (1975) 361. J.C. Myland and K.B. Oldham, Anal. Chem.. 60 (1988) 62. B.S. Grabaric, R.J. O’Halloran and D.E. Smith, Anal. Chim. Acta, 133 (1981) 349. S.O. Engblom and A.U. Ivaska in M.R. Smyth and J.G. Vos (Eds.), Analytical Chemistry Symposia Series, Vol. 25: Electrochemistry, Sensors and Analysis, Elsevier, Amsterdam, 1986, p. 49. I A. Erdelyi (Ed.), Tables of Integral Transforms, Vol. 1, McGraw-Hill, New York, 1954. p. 30. 8 J.C. Myland, K.B. Oldham and Zhu Guoyi, Anal. Chem.. 60 (1988) 1610. 9 I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, New York, 1965, p. 341. 10 D.F. Elliott and K.R. Rao, Fast Transforms: Algorithms, Analyses, Applications. Academic Press, Orlando, 1982, p. 38.
394 11 Ref. 10, p. 50. 12 D.E. Smith in A.J. Bard (Ed.), Electroanalytical Chemistry, Vol. 1, Marcel Dekker, New York, 1966, Ch. 1. 13 S. Bochner, Lectures on Fourier Integrals, Princeton University Press, Princeton, 1959, pp. 9, 10. 14 Ref. 9, p. 344.
J. Electroanal. Chem., 2% (1990) 371-394 Elsevier Sequoia S.A., Lausanne
The Fourier transform of a voltammetric peak and its use in resolution enhancement * Sten 0. Engblom Laboratory of Analytical
Chemistry, A*bo Akademi,
Biskopsgatan
8, SF-20500 Abe (Finland)
(Received 18 May 1990)
Abstract
The theory presented shows how the Fourier transform of a voltammetric peak can be utilized for finding peak parameters and enhancing the resolution. It is based on the assumption that the voltammettic or polarographic peak can be described by the se&-function. Potential uses of the Fourier transform of this function are discussed, but the main theme is resolution enhancement. New deconvolution-convolution functions for this purpose are introduced. It is possible to retain either the height of the original peak, or the area under it, by choosing the correct deconvolution function.
INTRODUCTION
One of the reasons, besides improved sensitivity, for introducing polarographic techniques with peak-shaped responses was originally to enhance the resolution when several electroactive species were present [l]. The possibilities for a simultaneous determination of several components are limited in classical dc polarography by the absence of a reliable base-line. In ac polarography, and other methods with similar responses, this difficulty is in many cases reduced. Also, these techniques are better suited to handle a situation where the different half-wave potentials are close to each other. However, in some cases serious overlapping of peaks does occur and a simultaneous determination is difficult, if not impossible. Several examples of this situation exist, among them solutions containing lead and thallium, or cadmium and indium. Many efforts have been made to either prevent the overlapping by chemical methods, or to improve the recorded curve using mathematical means. The subject of this article is the manipulation of the recorded data in the frequency domain, i.e. after Fourier transformation. Especially, the possibility of changing the width of the peaks without affecting their area, height or position, is
l
Presented at the J. Heyrovsky Centennial Congress on Polarography, Prague, 20-25
0022-0728/90,‘$03.50
Q 1990 - Elsevier Sequoia S.A.
August 1990.
372
treated. This operation, usually a deconvolution, is one of the major applications of Fourier transformation in ac polarography, according to Smith [2]. The theory presented below applies to reversible ac polarography, but the same type of response occurs in several electroanalytical techniques, e.g. semidifferential electroanalysis [3] and derivative neopolarography [4]. Grabaric et al. [5] introduced this type of resolution enhancement in electroanalytical chemistry and their deconvolution function was later slightly modified by Engblom and Ivaska [6]. However, this early work lacked the proper theory needed to explain the operation in detail. It is the purpose of this work to give that theoretical background, and at the same time suggest some new deconvolution and also convolution functions of use in resolution enhancement. All calculations in this work were performed using the MATLAB program (trademark of The MathWorks Inc.) on a VAX 8800 computer (trademark of Digital Equipment Corporation). The same program was used to generate the figures shown. THE CONTINUOUS
The general
CASE
form of the function
p(x)=cosh2 [b;lx
under
consideration
here is
(1)
- c)]
The main characteristics of this peak-shaped curve, see Fig. lA, are the height a, the width-determining parameter b, and the position c. A usual measure of width is the width at half height of the peak, Wl,2 = ln(3 + \/8)/b. Further, it is easily shown that the area under the curve equals 2a/b. For later reference, please note that a peak with a large b is narrower and the area under the curve is smaller in comparison with a peak with a smaller value for b. The Fourier transform of this function is given by [7] P(y)
=jm
p(x)
exp(-iny)dx=exp(-icy)?
sitiy/tJ2b)
--m Three main qualities of the peak have thus been separated into three factors: The first factor exp( - icy) is a function of the peak’s original position c, the second is the area under the peak, and the third factor is a function of the width of the peak, this time slightly above the half height, W,.,,,,,,. = r/2b. In Fig. 1B transforms are shown for different b’s. Here one can recognize a fact, which one feels intuitively is correct, namely that a sharp peak (large b) has a “wide” transform. Fast changes in a time-domain curve correspond to high frequencies. On the other hand, a very wide peak has a “sharp” transform. Very few, low frequencies are needed to describe the transform of such a curve. The main reason for performing a Fourier transformation is generally to separate signal from noise. Usually, as in this case, the signal is composed of mainly low frequencies and high-frequency noise is well separated. From eqn. (2) it is clear that not only is the noise separated from the signal, but also that some interesting
313
X
1
I
2-
B 1.5 z
l-
0.5 -
,, ...’
0-10 ----.---8
_--6
-4
-2
0
-._ ---_____ -_____ 2
4
6
8
10
Y Fig. 1. Three voltammetric peaks of different widths (A) and their corresponding Fourier transforms (B). The peak parameters were as follows: o = 1 and c = 5 for all peaks, and 6 = 1 ( -),2(------)and 3 (. . . . . .).
374
characteristics of the peak are disentangled. If a system exhibits only one peak it would be a very simple task for, e.g., an automated instrument to find the area under the peak using P(0) = 2a/b, and the position of the peak using cy = n7~ + y # 0, n = . . . , arctan{ -Im[P(y)]/Re[P(y)]), (2n - l)a/2 d cy 6 (2n + l)n/2, Here Im[ ] and -2, -1, 0, 1, 2 )...) the n = 0 case being the most straightforward. Re[ ] denote the function’s imaginary and real parts, respectively. These parameters are usually the only ones required in a routine analysis. Note that a conventional measurement of peak height and position in many cases involves the use of only one point on the curve. A more elaborate method would use some kind of curve fitting in the region of interest. In the procedure described above no fitting is needed and all recorded data points are used. Furthermore, the data used would be essentially noise-free. This approach has been used previously in the analysis of voltammetric waves [S]. However, the main theme of this article is the manipulation of the third factor of the right-hand side in eqn. (2). This factor determines the width of the time-domain peak and is readily accessible in the frequency domain. Multiplication with a suitable function would clearly be the most straightforward way to accomplish a desired change in width. It would thus be possible to affect the width without any knowledge of either a or c. In earlier attempts to change the width of the peak in this way, b was also assumed unknown. This, however, produced some unwanted side-effects. These early methods will be reviewed to expose the problem. No knowledge of any peak parameters available The multiplication of two functions in the frequency domain is equivalent to a convolution operation in the time domain. Thus, if U(y), V(y) and W(y) are the Fourier transforms of U(X), v(x) and w(x), then the convolution U(X) = /@‘,v( z)w(x - z)dz and the inverse transform of U(y) = V( y)W( y) yield the same result. Finding v(x) when U(X) and w(x) are known is called a deconvolution. By necessity (simply because the transform, eqn. (2), was not known), the earlier methods were deconvolutions. Deconvolution functions were defined in the time domain and their transforms were calculated numerically. The transform of the voltammetric peak was divided by the transform of the deconvolution function and, if the form of the deconvolution function had been chosen suitably, a narrower peak was the result of an inverse transform. However, the height of the new peak was not predictable, it was a function of the width of the deconvolution function but the exact relation was not known. The deconvolution functions used were of the same functional form as the peak but the peak height was either arbitrary [5] or equal to N, the number of points recorded [6]. Here, as an example, let us consider a deconvolution function with a peak height equal to b’/2, where b’ is the width-determining parameter in the deconvolution function. We have thus dpa(x)
=
b’/2 cash* ( b’x )
the transform
of which
is given
by eqn. (2) by letting
c = 0 and
a = b’/2,
and
375
replacing b with b’. The name of the function is derived from deconvolution (‘d ‘) and voltammetric peak-shape (‘p’), while the superscript ‘a’ indicates that this particular function does not affect the area under the peak it is operating on. The reason for this is that the area under the curve defined in eqn. (3) is equal to one. Also, this deconvolution function will not affect the position of the peak, because it is centered over the origin. The transform of d@(x) consists, as P(y), of three factors, but now the first two are both equal to one. Capital letters are used to indicate Fourier transforms, so that, e.g., the transform of &P(x) will be denoted gPa( y). The deconvolution is accomplished by, in the frequency domain, forming the quotient 2a b’ sinh(ry/2b’)
RP”(y) =p(Y)/Dp”(Y) =exd-icy)~
b
sinh(1TY/2b)
Only the third factor has been affected by the operation, the first and second factors, determining the position of and area under the peak, are still intact. The inverse transform of the result RP(y) is given by [7]
rpaCx) = ‘iii’
cosh[2b(
sin( oh/b’) x - c)] + cos( ah/b’)
assuming b’ > b. ‘r ’ in the function name refers to “result”. M’=(y) lacks an inverse transform in the case 6’ < b (see the Appendix). The result of a deconvolution using a deconvolution function of the same shape as the voltammetric peak is thus also a peak. In this case the area under it is still the same as under the original peak, 2a/b, but the peak height is now a complicated function of b’, given by a tan(rb/2b’)/(ab/2b’). The width of the peak at half height is also, as could be it to equal (2/b)ln[cos(ab/2b’) expected, a function of b’; one finds + cos’( rb/2b’) + 11. The useful range of b’ has already been limited to b’ > b. It is now of interest to investigate the behaviour of rp”(x) for different values of b’ and especially the limits b’ 4 b and b’ + 00. As b’ + b it is evident that the width of the peak tends to zero while the height tends to infinity, while the area remains constant. Thus, we have rp”(x) -+ (2a/b)6(x - c ) as b’ + b, where 6(x) is Dirac’s delta function. The other limit is likewise interesting as one can show that rp”(x) --, p(x) as b’ + 00. The resulting peak will thus vary between a delta function, which must be seen as the ideal response for a resolution-enhancing operation, and the original peak itself, as b’ varies in the interval [b, cc). Of course, once the analytical form of P(y) is known, the possibility of using convolution instead of deconvolution, is opened up. In the deconvolution operation the transform of the deconvolution function appears in the denominator in the quotient shown in eqn. (4). Since DPa(y) is itself peak-shaped and decreases rapidly towards zero as 1y 1 increases, there can be difficulties in implementing this operation in practice. A convolution, i.e. a multiplication, using c-P”(y)
=
sinh( ry/2b’) vy/2b’
376
X
The two component peaks are also shown Fig. 2. A curve with two similar, overlapping peaks ( -). (- - - - - -). The peak parameters were: a = 1 and b = 1 for both peaks, and c = 4 for the first peak and c = 6 for the second.
would on the other hand decrease the risk for the noise-enhancement inherent in the division. However, since CP”(y) + cc as y + cc it is clear that this operation will also enhance noise if carried out without proper precautions. More will be said on these problems in a later section. There is no time-domain equivalent to CP(y) (see the Appendix). A system consisting of several peaks, see e.g. Fig. 2, does not present a problem if one assumes that the different responses are additive. That is, if a curve consisting of several peaks can be written as
mp(x) =Plb> +p*(x) + ... +P”(x> where each pi(x), i = 1,2,. . ., n, is given by eqn. (1) after inserting parameters cli, bi and ci. The linear property of the Fourier transformation ensures that the transform of mp( x) is given by
where each P,(y) is given by eqn. (2) by inserting the proper parameters. Now, a multiplication with U”(y) or a division with DP(y) will, simply because of the distributive nature of those operations, affect each peak individually. In the inverse transform the effect on the peaks is determined by their original width, in accordance with the relations given above (for an example of a situation with two similar
311
1.8,
X
Fig. 3. Theoretical results of a resolution enhancement of the curve shown in Fig. 2. The use of eqn. (3) will preserve the area under the peak (A) while eqn. (9) leaves the height unchanged (B). The curves shown are: The original curve (), the result using b’ = 2 (- - - - - -) and the result using b’ = 1.5 ( . . . . ..).
378
peaks see Fig. 3A). If complete separation which is unchanged, is a useful measure.
is accomplished
the area under
the peak,
The width-determining parameter b assumed known In reversible voltammetry knowing b corresponds to knowing n, the number of electrons participating in the electrode reaction. Note that in a one-peak system b is relatively easy to find using the third factor in eqn. (2). One could, e.g., use the relation x = In (x/y + 7 1 + (x/y) ), which is valid for y = x/sinh(x), and find x in an iterative manner. This way of determining b would benefit from the advantages mentioned earlier, all points on a curve would contribute and the influence of noise would be eliminated to a large extent. Once b is known, it is a simple matter to introduce a factor into the convolution-deconvolution functions to achieve constant peak height in the result. In this case, dp becomes dph(x)
=
Similarly,
CP”( y)
tan( nb/2b’) ?rb/2b’
b’/2
the convolution =
mb/2b’ tan( nb/2b’)
(9)
cash* (b’x) function
is given by
sinh( vy/2b’) vy,‘2 b’
The superscript ‘h’ indicates b’. The resulting peak
(10)
here that the peak height of the result is independent
sin( nb/b’) rph(x)
= tan(ri/2b’)
cosh[2b(x
- c)] + cos(rb/b’)
of
(11)
is, of course, of the same shape as rp”(x) but the peak height is now equal to a for all b’ while the area under the peak is a strong function of b’; an integration yields the result ma/b’ tan(ab/2b’). The width of the peak at half height is still the same as that of rp”(x), or (2/b) In [cos(sb/2b’) + cos*(ab/2b’) + 11. rph(x) behaves in the same way as rp”(x) as b’ becomes very large, i.e. we have that rph(x) + p(x) as b’ + co. The situation is different at the other limit. As b’ -B b we see that both the width of the peak as well as the area under it tend to zero, but the peak height remains constant. This could lead one to expect a limit consisting of a finite delta function at x = c. That this limit in fact cannot be reached is evident in the frequency domain where we have RPh(y)
= exp(-icy)
b’ sinh( Ty/2 b' ) b’ tan(yb/2b’>
b sinh( ry/2b)
(12)
and it is clear that RPh(y) + 0 for all y as b’ + b. A multi-peak system is more problematic in this case. In Fig. 3B the simple case of two similar peaks is used, but, unfortunately, it is highly unlikely that all peaks have the same width. One then has to choose one b and accept that only peaks of
319
this width will retain their true heights. As an example, consider a system consisting of n peaks. Each peak p,(x) is characterized by the parameters ai, b, and ci, i=12 7 ,.*., n. If b, is chosen to be used in the resolution-enhancing operation, we get up(x)
=
tan( nbj/2b’)
b’/2
nb,/2b’
cosh’( b'x )
(13)
or CP”( y) =
rbj/2 b’ tan( sbj/2b’)
If the effect on peak p;(x) r#*(X)
=
sinh( vry/26’) I ry/2b i #j, one arrives at the result
is studied,
abj b, tan( abj/2b’)
(14)
sin( n-b/b’) (15)
- c)] + cos(rb;/b’)
cosh[2b,(x
b, tan( mbJ2b’) = bi tan(rbj/2b’) The factor relating
(16)
Vh(4
rph*(x) and rph(x) /> 1
if bi > b,
= 1
if bi = b,
bj tan( ab/2b’)
< 1
if bi < b,
bi tm( rbj/2b’)
’ -B I --+O +oo
in eqn. (16) has the following
as b’ + 00 as b’-+b,, b,
properties
(17)
bi>b,
In the last two cases it is assumed that it is meaningful to let b’ approach either bi or b,. In a system consisting of several peaks b’ should always be greater than, or equal to, the largest of the bi’s. This is required in order to ensure that inverse transforms of all the convolved peaks really exist (see the Appendix). Thus, in a multi-peak system, deconvolution using dph(x) or convolution using CPh( y) will produce the correct result only if the b in the deconvolution-convolution function is the same as in the original peak. The resulting peak will be smaller if the b used is larger than the b of the original peak, and larger if b in the deconvolution-convolution function is smaller than that of the peak. A natural question at this stage is if there are any other possible deconvolutionconvolution functions. Consider eqn. (2) again. It is evident that if b is known it is a simple matter to replace it with b’ while preserving the original shape of the peak. That is, as before it is possible to change the width of the peak without affecting, or needing to know, the area under it (alternatively, the height) or its position. To this one can now add that it is also possible to preserve the functional form of the peak.
380
All one has to do is to multiply b sinh( ry/2
the transform
of the peak, eqn. (2) with the function
b)
08)
cR”( ‘) = b’ sinh( ?ry/2b’) or, equivalently, divide the transform with DRa( y) = l/CP( y). function does not have an inverse transform (see the Appendix). transform of DRa( y) is given by [7] sin( Irb/b’) b’ dracx) = G cosh(2bx) + cos( mb/b’)
The CZV’( y) The inverse
(19)
What we have here is nothing but the result of a deconvolution using u”“(x) shown in eqn. (5), but for a constant factor. Hence the naming of this function, the first convolution letter in the function name, ‘d ‘, ‘c’ or ‘r’, still refer to deconvolution, and result, respectively, while the second letter, ‘r ‘, indicates that this function actually is of the same form as the result of a deconvolution-convolution of “the first kind”. The result of the operation is here
RR”(y)
= P( y)CP(
This function
y) =
has an inverse
P( y)/DP( similar
y) = exp(
--icy)?sinl(/tJib,) (20)
to eqn. (l),
(21) This new peak is of the same shape as the original one; they have the same positions but the new peak is narrower by the factor b/b’. Furthermore, the height of the new peak is increased by the factor b’/b. The area under the peak, however, remains the same (hence the superscript ‘a’). Clearly, it is also possible to derive functions of “the second kind” that will preserve the height of the original peak. These are simply defined by CRh( y) = (b/b’)CF( y) and DRh( y) = l/CRh( y). If P(y) in eqn. (2) is now operated on by either of these two, the result is
RRh(y)
= exp( -icy)%
sitiy/t!..b,)
with the time-domain
equivalent
rrh(X) = cosh2[b;x
- c>j
(22)
(23)
which gives a peak of the same height as the original one, but narrower by the factor b/b’. The area under the peak is also decreased by the same factor. In an operation of “the second kind” b’ is still restricted to the interval [b, 00). A b’ smaller than b would not be forbidden but it would produce a peak wider than the original. Also, in this case the convolution function would have an inverse
381
transform while the deconvolution function would lack one. However, the normal case is a b’ which is larger than b; in fact the larger b’ the narrower will the resulting peak be. As b’ + co we see that rP(x) + (2a/b)6(x - c). The fact that RRh(y) tends to zero for all y in this limit shows that the finite delta function, which is taking the value a for x = c but is zero everywhere else, is not attainable in the time domain. Of course, in both cases the resulting peaks tend to p(x) as
b’ -+ b. In a system exhibiting several peaks, the case where b in the deconvolution-convolution function is incorrect has to be considered once again. In an operation of the second kind this applies to both types of functions, i.e. both the area- and the height-preserving function. Assume that the width-determining parameter of peak j is used while the effect on peak i is studied, i #j. If, as an example, CRa( y) is used, one obtains
2a
RR”*(y) =exp(-icy)7
rY/2bi sinh( ry/2bi)
bj sinh( rry/2 bj ) 6’ sinh( ry/2b’)
bj sinh( Ty/2 bj ) = bi sinh( lry/2b,)
(24)
(25)
RR”( ‘)
The case is even more complicated here than in eqn. (16) since the proportionality factor is a function of y. The main properties are ‘>l
bj sinh( vy/2bj) bi sWry/2bi)
if 6, > b,
= 1
if bi = bj
< 1
if bi < bj
’ + 1
(26)
asi yl -+O asly
+m,
b,
--,CCI aslyl
*co,
bi>bj
-+O
If bi -c b, there is no problem regarding the existence of an inverse transform; we have I RR”*(y) I G I RRa( y) I for all y. This is sufficient since absolute integrability ensures the existence of the inverse (see the Appendix) and /‘?, I RRa*( y) Id y G /?, 1RF(y) ldy = 4nab’/b [9]. However, so far no explicit expression for the time-domain equivalent of RR”*(y) has been found. Something can still be said though about V”*(X). It is the convolution of rP(x) (eqn. 21) and a function given by eqn. (19) if b and 6’ in that equation are replaced by b, and b,, respectively. These two functions are symmetrical and peak-shaped and their convolution yields a symmetrical peak. This peak is centered at x = c as the original peak and the area under it is still 2a/b,. The width of this peak will be determined by not only b’ but also 6,. If b, is very close to bi so that the proportionality factor in eqn. (25) is only slowly decreasing towards zero, it will be possible to get a relatively narrow peak by using a large b’. On the other hand, if the difference in 6, and b, is large, even a very large b’ will produce only a minor effect on the width of the resulting peak.
382
The proportionality factor is thus in a way working like a digital filter where one can affect the cut-off frequency by varying bj. Similar conclusions are valid when considering the use of a height-preserving function. In this case the relation between RRh*(y) and RRh(y) is given by b,’ sinh( ry/2 RRh*( y) =
b, )
b: sinh( ny/2 b, )
The properties of the proportionality eqn. (26), but still different enough of them. We have b,? sinh( ry/2bj)
(27)
RRh(y)
factor in eqn. (27) are similar to those given in to warrant a separate exposition of at least some
> b,/b,
if b, > bj
< b,/b,
if 6, < bj
(28)
b,? sinh( ry/2b,) i -b,/b,
aslyl
-0
The missing ones (the b, = bj case and the limits as ) y 1 -+ co) are the same as in eqn. (26). In the b, < b, case the existence of the inverse can, as above, be drawn from the fact that /‘E’, 1RRh*( y) Id y d (b/b,)/‘?, I RRh( y) Id y = 4nab,/b, [9]. The shape of rrh’(x) is similar to that of rr”/*(x), but the two functions are scaled differently; we have rr”*(x) = (bb,/b’bi)rra*(x). If b, > 6, we have to look more closely at the behaviour of the frequency-domain function. As has been stated several times, absolute integrability guarantees the existence of an inverse. If the function under investigation is even (as RR” *( y ) and RRh*( y) are) the consideration can be restricted to the interval [0, cc). Furthermore, the integral can be divided into two parts, lo” ] g(x) I dx = 10”] g(x) I dx + /p ] g(x) I dx. If g(x) in this example is continuous on the closed interval [O,k] then ] g(x) ] is integrable. It is then the second integral which is of interest. When either RRa*( y) or RRh*( y) is considered, assume a k which is large enough to justify the approximation sinh( x) = exp(x). The integrand is then simplified and its behaviour for large y is given by (the constant factors are dropped and c = 0 is assumed for simplicity)
----7
= ”
i
0
if 6, -c bj
0
if b’ < b,b,/( bj - 6,)
and
bi > bj
00
if b’>b,bj/(bi-b,)
and
b,>b, (29)
It is clear that an integral will converge in the first two cases but not in the last case. If b, > b, the magnitude of 6’ is severly limited, which is a considerable disadvantage when considering that in a resolution-enhancing operation of the second kind a large 6’ means better resolution. In practise, dealing with experimental data one is of course obliged to use a discrete version of the results presented above. This will be the theme of the next section.
383
THE DISCRETE
CASE
Assume that an experimental are given by Pk = cosh2[b(xi:+
kh - c)]
curve is known
k=0,1,2
in N equidistant
points.
N-l
)...)
The points
(30)
where xin is the position of the first measured point and h is the distance between consecutive points. These points are of course related to the continuous function of the last section through pk =p(xi,, + kh). It is assumed here and below that the number of points, N, is large enough to characterize the peak and its Fourier spectrum completely. A specific number is discussed below. The discrete Fourier transform is defined by N-l
PI = c
pk exp( - i2nkl/N)
i=O,
1, 2 ,...,
N-
1
(31)
k=O
and the original
data are recovered
through
the inverse
transform
N-l
1 pk = x c P, exp(i2rkl/N) I=0
k=O,
1,2 ,...,
N-l
(32)
The discrete Fourier transform of the array of data in eqn. (30) should related to P(y) in eqn. (2). A reasonable model is P, = exp( - i2rQ/N)
7~‘l/bhN sinh( r’l/bhN
)
I=O,
1,2 ,...,
[N/2]
be closely
(33)
In the discrete case the position of the peak is taken relative to the first sample, i.e. if we have c = xi,, + .$h, 0 G [ G N - 1, then the first factor is a function of E instead of c. The second factor is now the sum of the recorded points rather than the area under the curve, while the third factor is the same as in eqn. (2) if y is replaced by 2rl/hN. It is known that, since pk is a real array, PN_, = P,* (the complex conjugate of P,) [lo]. The points given in eqn. (33) are then enough to generate the whole array. In eqn. (33) [N/2] signifies the greatest integer less than or equal to N/2, which is N/2 if N is even and (N - 1)/2 if N is odd. When a numerically calculated transform is compared with eqn. (33) the fit is perfect except for the highest frequencies, see Fig. 4. The reason for this deviation is not yet known. Equation (33) can now be used to estimate the minimum number of samples required. The third factor in eqn. (33) is rapidly decreasing towards zero. Only P, of small values on I are needed to describe the peak in the Fourier domain. Assume that PI with I> A does not contribute significantly to the transform, and that 1P,, 1 ( O.OOOlP,. The problem is then to find those values of z that satisfy z/sinh(z) < 0.0001. The answer is z >, 12.4231. In eqn. (33) z is r21/bhN. The product hN can be seen as the interval in x used when sampling the function and let us assume this is equal to five peak-widths, or hN = 5 ln(3 + 6)/b. This gives
384
-21 0
20
40
60
80
100
120
20
40
60
80
100
120
-l6ot.
v Fig. plot vs. h=
and eqn. (33) (....e.). In a normal 4. Comparison of a numerically calculated transform ( -) of 1P, 1 vs. I the fit seems perfect (A), while a small discrepancy shows up if, e.g., ln( 1P, I) is plotted I (B). The peak parameters were: a = 1, b = 3 and c = 5. The peak was sampled using xin = 0, 0.078125 and N = 128. The peak is shown in Fig. 1 (. . . . . .).
385
I = 62.1155 ln(3 + a)/;* = 11.0941. The same number of data is needed for the above-N/2 part of the array, so this crude analysis gives a lower limit for the number of samples at about N = 25. Of course, in order to sample high-frequency noise correctly in accordance with the sampling theorem, a much higher sampling rate is usually required. As in the continuous case, it is now possible to find E, and thus c, by using the relation 2a&/N = n?~ + arctan { -Im[P,]/Re[ P,]} which is valid when (2n 1) N/4 Q 51~ (2n + 1) N/4, I f 0, n = . . . , - 2, - 1, 0, 1, 2,. . . . Once again, n = 0 simplifies the calculation even though it restricts the magnitude of (1 considerably. Actually, even if only PI is used, 5 is still restricted to the first quarter of N. This restriction can be eased by adding zeros to the end of the data array before Fourier transformation. PO is useful when an estimate of the area under the peak is required. We have that P,, = If:: pk and therefore hP,, = 2a/b. Here it is necessary that the peak is entirely within the sampled region (meaning that the peak’s contribution to samples outside this region is negligible). It is also possible to find b as explained earlier. Resolution enhancement is achieved in the same manner and using essentially the same functions as in the previous section. A convolution is, in the discrete case, defined as a sum, uk = Cf:i u, w,_,, and the equivalent operation in the frequency domain is a term-by-term multiplication of the discrete Fourier transforms, U, = YW, [ll]. A deconvolution is then defined as a term-by-term division of two Fourier transforms. If a time-domain deconvolution function is used, one has to consider the periodicity implicit in eqns. (31) and (32). Both pk and P, are seen as periodic and, in order to avoid this periodic repetition being “seen” in a arrays, convolution-deconvolution operation, one has to destroy the periodicity. This is accomplished by zero-filling. An aperiodic array is created by appending N zeros to the original array [ll]. In eqn. (3) a deconvolution function was given that did not affect the area under the peak nor was any knowledge of the peak parameters required. In the discrete case this function has to be slightly modified in order to produce the same end result. The discrete transform is, as the continuous one, a product of three factors, see eqn. (33). In the transform of the deconvolution function it is necessary for the first two factors to equal one. The first factor has to equal one, i.e. 5 = 0, in order not to affect the original peak position. The maximum of the deconvolution function has then to occur at x = xi,, instead of at x = 0 as in the continuous case. The second factor, the sum of all points, is related to the area under the curve. If this sum is multiplied by h one gets the Riemann sum, which, if h is small enough, gives a good approximation to the area under the curve. For the deconvolution function not to affect the sum, and thereby the area under the peak, its own sum has to equal one. We can now define dp,” =
hb’/2 K=[-;+I]
,...,
-l,o,l,...,
I;]
cosh*( b’Kh)
and by using
dpY, = dpi_,
(this can be deduced
from eqn. 32) one can construct
386
must be done to the the complete array d’i, k = 0, 1, 2,. . . , N - 1. Zero-filling array in eqn. (34) before rearrangement. The relation to the continuous function is dp,” = h dp”(x, + oh). The area under the curve h dp”(x) is equal to h, and therefore, approximating the area under dpz with hC, dpi, it follows that C, dp,” = 1. The Fourier transform of dpl is then given by DP,a =
r21/b’hN sinh( r’l/b’hN)
1=0,1,2
,...,
[N/2]
(35)
with the reservation, as mentioned earlier, that the fit at the highest frequencies is not perfect. The result of a deconvolution using dpt yields a result similar to that in the continuous case. We get b’ sinh( a’l/b’hN)
RPP = exp( - i2r@/N)
b sinh( a21/bhN )
l=O,1,2
,...,
[N/2] (36)
or after an inverse 2ab’ rpi = rb cosh[26(
transformation sin( rb/b’) xi” + kh - c)] + cos( ah/b’)
k=O,
1, 2 ,...,
N-l
As in the continuous case, b’ can now vary in the interval [b, co). Equation describes rpt accurately for intermediate values of b’. However, as b’ tends there are deviations from this equation. In the limit we have
(37) (37) to b
if k = .$ (38) otherwise where 6,, is Kronecker’s delta function. In the other limit, b’ + 00, we have as before that rpt --, pk. In Fig. 5 the height of and area under rpi are compared with the height and area of the result in a numerical calculation. The fit is quite good, only a minor difference in peak height for 6’ close to b being detected. The peak height in this region is actually growing a little faster than predicted, before attaining a finite value at b’ = b. Note that a large b’ gives a very sharp dp,“, thus necessitating a large number of samples (small h) in order to achieve a correct transform. This would not be necessary if the operation started off directly with DPF or CP: = l/DPp. Using the present theory this is of course possible only if the highest frequencies do not contribute to the signal in any significant way. The other relevant properties of the peak in eqn. (37) were given in the previous section. It is now possible to examine a little more closely the effect of an early deconvolution function, N/cosh’(b’Kh) [6]. The peak height of this function should, however, be discussed first. In the context where this function was defined, slightly different versions of eqns. (31) and (32) were used. The discrete Fourier transform was defined as P, = (l/N& pk exp( -i2nkl/N) while the inverse was calculated according to pk = C, P, exp(i2nkl/N). The difference lies in the position-
ing of the factor l/N. This might be thought of as an insignificant distinction, but it is of some importance. Let once again U,, V, and W, be the discrete transforms of uk, vk and wk. The equivalent operation to U, = V,W, is then uk = (l/N)C, v,w,_,. If, however, eqns. (31) and (32) are used, it is uk = C, v,w,_,. The presence or absence of the factor l/N determines the peak height of the deconvolution function. Therefore, using the conventions of this article, this early deconvolution function is given by dp; =
1
k=[-;+l]
cosh’( b’kh)
The transform
)...,
)...,
is given by (with the same reservation m21/b’hN sinh( 1~‘l/b’hN
Forming
-1,&I
P,/DP,’
the quotient
RP,’ = exp( - i2n[f/N
f=O,
)
(39)
as earlier) IN/21
(40)
gives
) “;
) (crp;
1,2 ,...,
I;]
sffL[ ~~~~~~)’
I = 0, 1, 2,. . . ,[ N/2]
(41)
k
where Crp: = Cp,/Edpi. rp; =
given by
sin( sb/b’)
a( b’)2h mb
The inverse
cosh[2b(
xi” + kh - c)] + cos( ah/b’)
k=O,
1,2 ,...,
N-l (42)
again applies only for intermediate values of b’. As b’ -+ b we have that rpz + as,*, instead of the limiting infinite delta function indicated by eqn. (42). At the other limit (b’ + 00) we have, for a fixed N (and h), that dp: + Sk,, the transform of which equals unity for all 1. A deconvolution has no effect and the result equals the original. In the case where b is known to the experimenter, it is also, in the discrete case, possible to use deconvolution-convolution functions that will preserve the peak height and/or the functional form of the original peak. The discrete version of eqn. (9) is given by dp, =
tan( rb/2b’) nb/2b’
dp” K
.=I-;+1]
,...,
As before, the range of b’ where this function the use of this function is rpkh= tan(rij26’) k=O,
1,2 ,...,
cosh[2b(xi, N-l
-l,o,l,...
:I;]
can be used is limited.
The result of
sin( nb/b’) + kh - c)] + cos(rb/b’) (44)
18 -
A
16 14 12 . 10 86-’ 42-
* . ‘* l**.
2
l *.*.
3
l
4
. 5
.
*
.
l
.
6
1
8
9
10
“--‘-1.2 -
B
I*- . . . . . .t.
. . . . .
l
.
.
*
l
. i
0.8 -
0.6 -
0.4 -
0.2 -
2
3
4
5
6
7
8
9
10
b’
Fig. 5. Resulting peak height (A) and area under the peak (B) as a function of the b’ used in a deconvolution. The predicted height and area of the result rpi (eqn. 37) (. . . . . .) is compared with the result in an actual numerical calculation ( l). In the numerical example eqn. (30) was used to generate synthetic data. The peak parameters were: (I = 1, 6 = 2, c = 5, xin = 0, h = 0.078125 and N = 128 (256 after zero-filling). Equation (34) was used to construct the deconvolution function. The simple formula hE, rp; wasusedto calculate the areas after deconvolution.
389
1.4
.
A
1.24 l
1-
s
0.8 -
: &
0.6 -
..**.*. .
l
l
. .
.
.
l
.
l
.
l
8
9
0.4 -
0.2 -
01 2
3
4
5
6
7
b’
1.4. B
1.2 c
l. 0.8 -
.
*
.
*
.
l
I
8
9
l
I 0.61
l
;
I .
l
l
0.4 -
.
’
t 0.2 - ’ l
01 2
3
4
5
6
10
b’
Fig. 6. A comparison of peak height (A) and area under the peak (B) with the predicted values (. . . . . .). The height-preserving deconvolution function dp,h in eqn. (43) was used in a numerical calculation and the results ( l) were compared with e-qn. (44). The peak parameters were as in Fig. 5. The areas were calculated using hE, rp*“.
390
1.8 1.6
A
1.4
0.8
I-
O-
I
2
3
4
5
6
I
8
9
10
6
7
8
9
10
were generated according Fig. 7. Results of a deconvolution using drxa(eqn. 45). Synthetic data ( -) to eqn. (30) using D = 1, h = 3, c = 5, h = 0.078125 and N = 128 (256 after zero-filling). Results using b’ = 4 (- - - - - -) and b’ = 5 (. . . . .) are shown (A). A logarithmic analysis (B) shows the correct slopes long linear parts are due to the fact that (3, 4 and 5) and the correct intercept (x = 5). The exceptionally the original peak was computer-generated and noise-free.
391
See Fig. 6 for an example of how well the height of and area under the resulting peak in a real calculation adheres to the height and area of the peak in eqn. (44). The fit is again very good except for a deviation in height when b’ is very close to b. The height is, as before, growing faster until the tangent-factor in the denominator forces the peak height towards zero. Considering a multi-peak system, the discussion in the provision section still applies wherever eqn. (44) is valid. Discrete versions of the functions that preserve the functional form of the peak, are also available. They will be given as deconvolution functions below, but the convolution functions can also be used if one takes into account the deviating behaviour at the highest frequencies (e.g. by using filters). The area-preserving deconvolution function of the second kind is thus given by I drKa= c
sin( vrb/b’) cosh[ 2b( xin + id)]
+ cos( ah/b’)
.=[-;+I] ,...,-1,0,1,...,
[;J
(45)
1.8 peak 1
2
1.61.41.2 I-
“2 c
0.8 0.6 0.4 0.2 0
\t.. 0
1
2
3
4
5
6
I
8
9
I 10
X
Fig. 8. Deconvolution of a double-peak system. The peak parameters were for peak 1: a = 1, b = 2 and c = 4.5, and for peak 2: CI= 1, b = 3, and c = 5.5. The peaks were generated separately using eqn. (30) with h = 0.078125 and N = 128 (256 after zero-filling), and then added. In the deconvolution function b of peak 2 was used, i.e. b = 3. The results of using b’ = 4 (- - - - - -) and b’ = 5 (. . . . . .) are shown together with the original curve (). The resulting peak heights of peak 1 cannot yet be predicted. The peak heights of peak 2 are (the true values are shown within brackets): 1.071 (l.C!C@). 1.369 (1.333) and 1.674 (1.667).
392
and the result by ab’/b “r’=
cosh2[b’(xi,+kh-c)]
k=O,
1, 2 ,...,
N-l
The height-preserving function is simply related to dr,” through drf’ = (b’/b)dr:, and the result is similarly given by ori = (b/b’)rr,“. As before, eqn. (46) is valid for a wide range of b’. However, for very large b’ and for b’ close to b there are deviations. As 6’ + co we have that dr,” --+ (hb/2)/cosh2(b~h). We have thus for large 6’ that rrt is following closely the behaviour of rpi for b’ close to b. rrf tends to zero in this limit as a consenquence of the factor b/b’ and the fact that the limiting delta function is finite. In the less interesting limit, b’ + b, both rr: and rrt tend to pk according to eqn. (46). However, the use of eqn. (45) are not practical in this limit. In Fig. 7 some results of the use of eqn. (45) is shown. A logarithmic analysis [12] of the peaks was performed to show that the curves were of the right shape. Finally, Fig. 8 shows a double-peak system with different b’s and the response to a deconvolution using drKa. b of peak 2 was used in the deconvolution function and it is thus that peak which shows a predicted response to the deconvolution. The resulting peak heights of peak 1 are not yet predictable; they are however still directly proportional to a of that peak. The areas under both peaks remain the same and if complete separation is achieved, these can be used. CONCLUSIONS
In this work the theoretical background for resolution enhancement of voltammetric and polarographic peaks by deconvolution has been given. Utilizing the Fourier transform of a peak, several functions were derived, each with a distinct property distinguishing it from the others. The user can choose between two forms of the resulting peaks, and whether the original peak height or the area under the peak shall be retained. Numerical calculations using synthetic data yielded results that agreed well with the theoretical predictions. The present theory works well both in the continuous and the discrete case. Some refinement is still needed in the discrete case for the high-frequency region, but in practice, when filters are in any case needed to control instrumental and other noise, the theory given here is fully adequate. If compared with other methods of finding the components of a fused system of peaks, the present method has the advantage of not needing any initial guesses with regard to number of peaks, positions of peaks, peak heights, peak widths, etc. No iterations of any sort are needed in the resolution-enhancing operation, thus eliminating cases of lengthy, diverging calculations. Among the disadvantages, one can mention the fact that the present theory concerns only reversible peaks; the incorporation of quasi-reversibility into the theory can be very difficult. In the crucial deconvolution-convolution step a sensitivity to noise is characteristic.
393
APPENDIX
The inverse
transform
of a function
P(y)
is given by (Al)
Generally, the integral in eqn. (Al) exists for all x if P(y) is absolutely integrable, or at least for all x # 0, if P(y) converges monotonically to zero as y + f cc [13]. The functional form of interest here occurs in several of the definitive equations in this article. CRa( y) (eqn. 18) is chosen as an example, but the results below apply to all cases after obvious modifications. That this function is absolutely integrable when b’ < b is well known; the integral can be found in most handbooks, see e.g. ref. 14,
If b’ > b, the integrand above will grow beyond any limits as y + cc. To see this use, e.g., the approximation sinh(x) = exp(x) which applies for very large x. Thus, the function is absolutely integrable and an inverse will exist if and only if b’ < b. The inverse is in fact given by eqn. (19) by letting b and b’ change places. The reasoning above could equally well be applied to RPa( y ) or DRa( y ) and explain why the time-domain functions rp”(x) and d?(x) exist when b’ > b but do not when b’ -c b. The functions CRh( y), DRh( y) and RPh( y) differ from the above-mentioned functions only by a constant factor, and are therefore also subject to the same conclusions regarding the existence of their time-domain equivalents.
REFERENCES
1 G.C. Barker in P. Zuman (Ed.), Progress in Polarography. Vol. 2, Interscience, New York, 1962, Ch. 19. 2 D.E. Smith in A.G. Marshall (Ed.), Fourier, Hadamard, and Hilbert Transforms in Chemistry, Plenum Press, New York, 1982, Ch. 15. M. Goto and D. Ishii, J. Electroanal. Chem., 61 (1975) 361. J.C. Myland and K.B. Oldham, Anal. Chem.. 60 (1988) 62. B.S. Grabaric, R.J. O’Halloran and D.E. Smith, Anal. Chim. Acta, 133 (1981) 349. S.O. Engblom and A.U. Ivaska in M.R. Smyth and J.G. Vos (Eds.), Analytical Chemistry Symposia Series, Vol. 25: Electrochemistry, Sensors and Analysis, Elsevier, Amsterdam, 1986, p. 49. I A. Erdelyi (Ed.), Tables of Integral Transforms, Vol. 1, McGraw-Hill, New York, 1954. p. 30. 8 J.C. Myland, K.B. Oldham and Zhu Guoyi, Anal. Chem.. 60 (1988) 1610. 9 I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, New York, 1965, p. 341. 10 D.F. Elliott and K.R. Rao, Fast Transforms: Algorithms, Analyses, Applications. Academic Press, Orlando, 1982, p. 38.
394 11 Ref. 10, p. 50. 12 D.E. Smith in A.J. Bard (Ed.), Electroanalytical Chemistry, Vol. 1, Marcel Dekker, New York, 1966, Ch. 1. 13 S. Bochner, Lectures on Fourier Integrals, Princeton University Press, Princeton, 1959, pp. 9, 10. 14 Ref. 9, p. 344.