spec.fft: 1D/2D/nD (multivariate) spectrum of the Fourier transform
In spectral: Common Methods of Spectral Data Analysis
Description
Usage
Arguments
Details
Value
Missing Values
See Also
Examples
Description
This function calculates the Fourier spectrum and power spectral density
of a given data object. The dimension of the array can be of arbitary size
e. g. 3D or 4D.
Usage
1
Arguments
y
1D data vector, y coordinate of a 2D matrix, nD (even 2D) array
or object of class fft
x
x-coordinate of the data in y or z. If y is an array, x must be a named list x = list(x = ..., y = ...).
z
optional 2D matrix
center
logical vector, indicating which axis to center in frequency space
Details
The function returns an user friendly object, which contains as much frequency
vectors as ordinates of the array. spec.fft provides the
ability to center the spectrum along multiple axis. The amplitude output is already
normalized to the sample count and the frequencies are given in terms of
1/Δ x-units.
Value
An object of the type fft is returned. It contains the
spectrum A, with "reasonable" frequency vectors along each ordinate. psd represents
the standardized power spectral density, [0,1]. The false alarm probability (FAP)
p is given similar to the Lomb-Scargle method, see spec.lomb.
Missing Values
Given a regualar grid x_i = δ x \cdot i there might be missing values
marked with NA, which are treated by the function as 0's.
This "zero-padding" leads to a loss of signal energy being
roughly proportional to the number of missing values.
The correction factor is then (1 - Nna/N) as long as Nna / N < 0.2.
If the locations of missing values are randomly
distributed the implemented procedure workes quite robust. If correalted
gaps are present, the proposed correction is faulty and
scales wrong. This is because a convolution of the incomplete
sampling vector with the the signal takes place. An aliasing effect
takes place distorting the spectral content.
To be compatible with the underlying Fourier transform, the amplitudes
are not affected by this rescaling.
Only the power spectral density (PSD) is corrected in terms of the energy
content, which is experimental for the moment.
See Also
plot.fft
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
# 1D Example with two frequencies
#################################
x <- seq(0, 1, length.out = 1e3)
y <- sin(4 * 2 * pi * x) + 0.5 * sin(20 * 2 * pi * x)
FT <- spec.fft(y, x)
par(mfrow = c(2, 1))
plot(x, y, type = "l", main = "Signal")
plot(
FT,
ylab = "Amplitude",
xlab = "Frequency",
type = "l",
xlim = c(-30, 30),
main = "Spectrum"
)
summary(FT)
# 2D example with a propagating wave
####################################
x <- seq(0, 1, length.out = 50)
y <- seq(0, 1, length.out = 50)
# calculate the data
m <- matrix(0, length(x), length(y))
for (i in 1:length(x))
for (j in 1:length(y))
m[i, j] <- sin(4 * 2 * pi * x[i] + 10 * 2 * pi * y[j])
# calculate the spectrum
FT <- spec.fft(x = x, y = y, z = m)
# plot
par(mfrow = c(2, 1))
rasterImage2(x = x,
y = y,
z = m,
main = "Propagating Wave")
plot(
FT,
main = "2D Spectrum",
palette = "wb"
,
xlim = c(-20, 20),
ylim = c(-20, 20),
zlim = c(0, 0.51)
,
xlab = "fx",
ylab = "fy",
zlab = "A",
ndz = 3,
z.adj = c(0, 0.5)
,
z.cex = 1
)
summary(FT)
# 3D example with a propagating wave
####################################
# sampling vector
x <- list(x = seq(0,2,by = 0.1)[-1]
,y = seq(0,1, by = 0.1)[-1]
,z = seq(0,1, by = 0.1)[-1]
)
# initializing array
m <- array(data = 0,dim = sapply(x, length))
for(i in 1:length(x$x))
for(j in 1:length(x$y))
for(k in 1:length(x$z))
m[i,j,k] <- cos(2*pi*(1*x$x[i] + 2*x$y[j] + 2*x$z[k])) + sin(2*pi*(1.5*x$x[i]))^2
FT <- spec.fft(x = x, y = m, center = c(TRUE,TRUE,FALSE))
par(mfrow = c(2,2))
# plotting m = 0
rasterImage2( x = FT$fx
,y = FT$fy
,z = abs(FT$A[,,1])
,zlim = c(0,0.5)
,main="m = 0"
)
# plotting m = 1
rasterImage2( x = FT$fx
,y = FT$fy
,z = abs(FT$A[,,2])
,zlim = c(0,0.5)
,main="m = 1"
)
# plotting m = 2
rasterImage2( x = FT$fx
,y = FT$fy
,z = abs(FT$A[,,3])
,zlim = c(0,0.5)
,main="m = 2"
)
rasterImage2( x = FT$fx
,y = FT$fy
,z = abs(FT$A[,,4])
,zlim = c(0,0.5)
,main="m = 3"
)
summary(FT)
# calculating the derivative with the help of FFT
################################################
#
# Remember, a signal has to be band limited.
# !!! You must use a window function !!!
#
# preparing the data
x <- seq(-2, 2, length.out = 1e4)
dx <- mean(diff(x))
y <- win.tukey(x) * (-x ^ 3 + 3 * x)
# calcualting spectrum
FT <- spec.fft(y = y, center = TRUE)
# calculating the first derivative
FT$A <- FT$A * 2 * pi * 1i * FT$fx
# back transform
dm <- spec.fft(FT)
# plot
par(mfrow=c(1,1))
plot(
x,
c(0, diff(y) / dx),
type = "l",
col = "grey",
lty = 2,
ylim = c(-4, 3)
)
# add some points to the line for the numerical result
points(approx(x, Re(dm$y) / dx, n = 100))
# analytical result
curve(-3 * x ^ 2 + 3,
add = TRUE,
lty = 3,
n = length(x))
legend(
"topright",
c("analytic", "numeric", "spectral"),
title = "diff",
lty = c(3, 2, NA),
pch = c(NA, NA, 1),
col=c("black","grey","black")
)
title(expression(d / dx ~ (-x ^ 3 + 3 * x)))
Example output
Loading required package: rasterImage
Loading required package: plotrix
Loading required package: lattice
Length Class Mode
fx 1000 -none- numeric
A 1000 -none- complex
mode 1 -none- character
center 1 -none- logical
Length Class Mode
fx 50 -none- numeric
fy 50 -none- numeric
A 2500 -none- complex
mode 1 -none- character
center 2 -none- logical
Length Class Mode
fx 20 -none- numeric
fy 10 -none- numeric
fz 10 -none- numeric
A 2000 -none- complex
mode 1 -none- character
center 3 -none- logical
spectral documentation built on March 29, 2021, 5:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Missing Values See Also Examples
Description
This function calculates the Fourier spectrum and power spectral density of a given data object. The dimension of the array can be of arbitary size e. g. 3D or 4D.
Usage
1 |
Arguments
y |
1D data vector, y coordinate of a 2D matrix, nD (even 2D) array
or object of class |
x |
x-coordinate of the data in |
z |
optional 2D matrix |
center |
logical vector, indicating which axis to center in frequency space |
Details
The function returns an user friendly object, which contains as much frequency
vectors as ordinates of the array. spec.fft provides the
ability to center the spectrum along multiple axis. The amplitude output is already
normalized to the sample count and the frequencies are given in terms of
1/Δ x-units.
Value
An object of the type fft is returned. It contains the
spectrum A, with "reasonable" frequency vectors along each ordinate. psd represents
the standardized power spectral density, [0,1]. The false alarm probability (FAP)
p is given similar to the Lomb-Scargle method, see spec.lomb.
Missing Values
Given a regualar grid x_i = δ x \cdot i there might be missing values
marked with NA, which are treated by the function as 0's.
This "zero-padding" leads to a loss of signal energy being
roughly proportional to the number of missing values.
The correction factor is then (1 - Nna/N) as long as Nna / N < 0.2.
If the locations of missing values are randomly
distributed the implemented procedure workes quite robust. If correalted
gaps are present, the proposed correction is faulty and
scales wrong. This is because a convolution of the incomplete
sampling vector with the the signal takes place. An aliasing effect
takes place distorting the spectral content.
To be compatible with the underlying Fourier transform, the amplitudes are not affected by this rescaling. Only the power spectral density (PSD) is corrected in terms of the energy content, which is experimental for the moment.
See Also
plot.fft
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 | # 1D Example with two frequencies
#################################
x <- seq(0, 1, length.out = 1e3)
y <- sin(4 * 2 * pi * x) + 0.5 * sin(20 * 2 * pi * x)
FT <- spec.fft(y, x)
par(mfrow = c(2, 1))
plot(x, y, type = "l", main = "Signal")
plot(
FT,
ylab = "Amplitude",
xlab = "Frequency",
type = "l",
xlim = c(-30, 30),
main = "Spectrum"
)
summary(FT)
# 2D example with a propagating wave
####################################
x <- seq(0, 1, length.out = 50)
y <- seq(0, 1, length.out = 50)
# calculate the data
m <- matrix(0, length(x), length(y))
for (i in 1:length(x))
for (j in 1:length(y))
m[i, j] <- sin(4 * 2 * pi * x[i] + 10 * 2 * pi * y[j])
# calculate the spectrum
FT <- spec.fft(x = x, y = y, z = m)
# plot
par(mfrow = c(2, 1))
rasterImage2(x = x,
y = y,
z = m,
main = "Propagating Wave")
plot(
FT,
main = "2D Spectrum",
palette = "wb"
,
xlim = c(-20, 20),
ylim = c(-20, 20),
zlim = c(0, 0.51)
,
xlab = "fx",
ylab = "fy",
zlab = "A",
ndz = 3,
z.adj = c(0, 0.5)
,
z.cex = 1
)
summary(FT)
# 3D example with a propagating wave
####################################
# sampling vector
x <- list(x = seq(0,2,by = 0.1)[-1]
,y = seq(0,1, by = 0.1)[-1]
,z = seq(0,1, by = 0.1)[-1]
)
# initializing array
m <- array(data = 0,dim = sapply(x, length))
for(i in 1:length(x$x))
for(j in 1:length(x$y))
for(k in 1:length(x$z))
m[i,j,k] <- cos(2*pi*(1*x$x[i] + 2*x$y[j] + 2*x$z[k])) + sin(2*pi*(1.5*x$x[i]))^2
FT <- spec.fft(x = x, y = m, center = c(TRUE,TRUE,FALSE))
par(mfrow = c(2,2))
# plotting m = 0
rasterImage2( x = FT$fx
,y = FT$fy
,z = abs(FT$A[,,1])
,zlim = c(0,0.5)
,main="m = 0"
)
# plotting m = 1
rasterImage2( x = FT$fx
,y = FT$fy
,z = abs(FT$A[,,2])
,zlim = c(0,0.5)
,main="m = 1"
)
# plotting m = 2
rasterImage2( x = FT$fx
,y = FT$fy
,z = abs(FT$A[,,3])
,zlim = c(0,0.5)
,main="m = 2"
)
rasterImage2( x = FT$fx
,y = FT$fy
,z = abs(FT$A[,,4])
,zlim = c(0,0.5)
,main="m = 3"
)
summary(FT)
# calculating the derivative with the help of FFT
################################################
#
# Remember, a signal has to be band limited.
# !!! You must use a window function !!!
#
# preparing the data
x <- seq(-2, 2, length.out = 1e4)
dx <- mean(diff(x))
y <- win.tukey(x) * (-x ^ 3 + 3 * x)
# calcualting spectrum
FT <- spec.fft(y = y, center = TRUE)
# calculating the first derivative
FT$A <- FT$A * 2 * pi * 1i * FT$fx
# back transform
dm <- spec.fft(FT)
# plot
par(mfrow=c(1,1))
plot(
x,
c(0, diff(y) / dx),
type = "l",
col = "grey",
lty = 2,
ylim = c(-4, 3)
)
# add some points to the line for the numerical result
points(approx(x, Re(dm$y) / dx, n = 100))
# analytical result
curve(-3 * x ^ 2 + 3,
add = TRUE,
lty = 3,
n = length(x))
legend(
"topright",
c("analytic", "numeric", "spectral"),
title = "diff",
lty = c(3, 2, NA),
pch = c(NA, NA, 1),
col=c("black","grey","black")
)
title(expression(d / dx ~ (-x ^ 3 + 3 * x)))
|
Example output
Loading required package: rasterImage
Loading required package: plotrix
Loading required package: lattice
Length Class Mode
fx 1000 -none- numeric
A 1000 -none- complex
mode 1 -none- character
center 1 -none- logical
Length Class Mode
fx 50 -none- numeric
fy 50 -none- numeric
A 2500 -none- complex
mode 1 -none- character
center 2 -none- logical
Length Class Mode
fx 20 -none- numeric
fy 10 -none- numeric
fz 10 -none- numeric
A 2000 -none- complex
mode 1 -none- character
center 3 -none- logical
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.