Fourier¶

CalcpadCE treats both classical analytical Fourier series and the discrete fast Fourier transform as first-class operations, exposing every coefficient and spectrum on the report instead of inside opaque library calls.

The analytical expansion computes Fourier coefficients of a user-defined function via direct integration and reconstructs it for any chosen number of terms. The classic-vs-FFT comparison puts the same function side-by-side under both methods, making the equivalence and aliasing limits of FFT-based reconstruction visible. The synthetic-signal FFT sweeps a two-tone signal through fft() and recovers power and amplitude spectral densities to illustrate the basic vibration-analysis workflow. The strong-motion record applies the same pipeline to a real earthquake acceleration history (Anamizu, 2024) to extract its Fourier amplitude spectrum.

Sampling rate, signal length and retained-term count are parametric, so each worksheet doubles as a sandbox for convergence, leakage and resolution studies.

FFT¶

Discrete-Fourier analysis of a synthetic two-tone signal that switches frequency mid-record, using the built-in fft() function. Returns the power and amplitude spectral densities so the basic vibration-analysis workflow stays fully visible on the worksheet.

Code:

#rad
'Sampling frequency -'Fs = 1024Hz
'Number of sampling points - 'N = 1024
'Time step -'Δt = 1/Fs
'Times of sampling points -'t = range_hp(0; N - 1; 1)*Δt
'Frequencies and amplitudes
f_1 = 10Hz','f_2 = 30Hz
X_1 = 20mm','X_2 = 20mm
t_0 = 0.5*N*Δt
f(t) = f_1*(t < t_0) + f_2*(t ≥ t_0)
X(t) = X_1*(t < t_0) + X_2*(t ≥ t_0)
'Input signal -'x(t) = X(t)*sin(2*π*f(t)*t)
x = hp([x(t)|0])
'Solution -'f = fft(x)
'Squared magnitude -'S = row(f; 1)^2 + row(f; 2)^2
'Power spectral density -'PSD = S/(N*Fs)
'Amplitude spectral density -'ASD = sqrt(PSD)
'<!--'PlotWidth = 600','PlotHeight = 150','PlotAdaptive = 0'-->
'Input signal plot
$Plot{line(i; row(x; 1)) & spline(i; row(x; 2)) @ i = 1 : N}
'Power spectral density plot
$Plot{spline(i; PSD) @ i = 1 : N/2}

Fourier Series¶

Analytical Fourier-series expansion of a user-defined function on a finite interval, computing each coefficient by direct numerical integration. Reconstructs the function for any chosen number of terms and overlays the truncated series against the original to show convergence.

Code:

#rad
'Interval length -'l = 10
'Function -'f(x) = (e^x - x^2 - 2*x^4)/l
'Coefficients -'a_0 = 2/l*$Integral{f(x) @ x = 0 : l}
a(k) = 2/l*$Integral{f(x)*cos(2*π*k/l*x) @ x = 0 : l}
b(k) = 2/l*$Integral{f(x)*sin(2*π*k/l*x) @ x = 0 : l}
'Expnsion to series
f_s(x) = a_0/2 + $Sum{a(k)*cos(2*π*k/l*x) + b(k)*sin(2*π*k/l*x) @ k = 1 : n}
'Plot for'n = 20'iterations
$Plot{f(x) & f_s(x) @ x = 0 : l}

Fourier Series + FFT¶

Side-by-side comparison of the same function expanded by classical Fourier-coefficient integration and by the discrete FFT. Highlights where the two reconstructions agree and how the sample count \(n\) controls FFT-based aliasing.

Code:

"Fourier series expansion (classic)
#rad
'Number of iterations -'n = 20
'Interval length -'l = 10
'Function -'f(x) = (e^x - x^2 - 2*x^4)/l
'Coefficients -'a_0 = 2/l*$Integral{f(x) @ x = 0 : l}
a(k) = 2/l*$Integral{f(x)*cos(2*π*k/l*x) @ x = 0 : l}', ' _
b(k) = 2/l*$Integral{f(x)*sin(2*π*k/l*x) @ x = 0 : l}
'Expansion to series
f_s(x) = a_0/2 + $Sum{a(k)*cos(2*π*k/l*x) + b(k)*sin(2*π*k/l*x) @ k = 1 : n}
'Plot the function and Fourier approximation
$Plot{f(x) & f_s(x) @ x = 0 : l}
"Fourier series expansion (FFT)
'<hr/>
'Number of samples -'n = 32
'Sampling step -'dx = l/n
'Sample points -'x = range(0; n; 1)*dx
'Apply FFT -'y = fft(hp([f(x)|0]))
'Fourier coefficients from FFT:
a_0 = y.(1; 1)/n'- first coefficient
a(k) = y.(1; k + 1)/n'- cosine coefficients, ' _
b(k) = -(y.(2; k + 1)/n)'-sine coefficients
x_t(x) = l - x/2'- abscissa scalling
'Function reconstruction
f_fft(x) = a_0/2 + $Sum{a(k)*cos(2*π*k/l*x_t(x)) + b(k)*sin(2*π*k/l*x_t(x)) @ k = 1 : n}
'Plot the function and its FFT-based Fourier approximation
$Plot{f(x) & f_fft(x) @ x = 0 : l}

Strong Motion FFT¶

FFT analysis of the Anamizu 2024 earthquake strong-motion acceleration record. Extracts the Fourier amplitude spectrum and power spectral density from a real seismic time history to identify dominant frequency content.

Code:

'Accelleration data
'Raw record length = 225.000 sec, Uncor max = 136992, at 73.050 sec.
'Processed: 2024-01-05 23:22:02 GMT, USGS, Max = 7.283 cm/sec2 at   72.590 sec
'Record filtered below 0.10 Hz (periods over 10.00 secs), and above 40.0 Hz
#read a_g from Anamizu Earthquake Strong Motion.txt TYPE=R_hp SEP=';'
'Number of record points -'n = n_cols(a_g)
'Time function 't(k) = (k - 1)/100*s
'Input accelleration, cm/s² -'a_g
'Sampling frequency -'Fs = 100Hz
'Frequency step -'Δf = Fs/n
'Frequency function -'f(k) = (k - 1)*Δf
'Number of output points -'N = floor(n/2) + 1
'FFT analysis -'A_f = submatrix(fft(a_g); 1; 2; 1; N)
'Real values -'Re = row(A_f; 1)
'Imaginary values -'Im = row(A_f; 2)
'Squared magnitudes -'S = Re^2 + Im^2
'Fourier amplitude spectrum (FAS) -'FAS = sqrt(S)/n
'Power spectral density -'PSD = S/(Fs*n)
'Magnitude spectral density -'ASD = sqrt(PSD)
'<!--'PlotHeight = 120','PlotSVG = 1','PlotAdaptive = 0'-->
'Input accelleration plot [cm/s²] vs time [s]
$Plot{t(k)|spline(k; a_g) @ k = 1 : n}
'Fourier amplitude spectrum (FAS) [cm/s] vs frequency [Hz] in the range of 0-10 Hz
$Plot{f(k)|spline(k; FAS) @ k = 1 : N/5}

Rendered Output:

Accelleration data

Raw record length = 225.000 sec, Uncor max = 136992, at 73.050 sec.

Processed: 2024-01-05 23:22:02 GMT, USGS, Max = 7.283 cm/sec2 at 72.590 sec

Record filtered below 0.10 Hz (periods over 10.00 secs), and above 40.0 Hz

Matrix a_g was successfully read from Anamizu Earthquake Strong Motion.txt SEP=';' TYPE=R

Number of record points - n = n_cols ( a_g ) = 22500

Time function t ( k ) = k − 1100 · s

Input accelleration, cm/s² - a_g = 7.18×10^-61.01×10^-59.21×10^-69.73×10^-69.41×10^-69.67×10^-69.5×10^-69.65×10^-69.57×10^-69.65×10^-69.6×10^-69.66×10^-69.66×10^-69.73×10^-69.76×10^-69.82×10^-69.89×10^-69.94×10^-69.99×10^-610^-5⋯0.0284

Sampling frequency - Fs = 100 Hz

Frequency step - Δf = Fsn = 100 Hz22500 = 0.00444 Hz

Frequency function - f ( k ) = ( k − 1 ) · Δf

Number of output points - N = floor(n2) + 1 = floor(225002) + 1 = 11251

FFT analysis - A_f = submatrix ( fft ( a_g ) ; 1; 2; 1; N ) = submatrix ( fft ( a_g ) ; 1; 2; 1; 11251 ) = 1.06-0.279-0.9420.8090.555-1.210.1121.29-0.952-0.8511.65-0.148-1.771.460.938-2.751.034.09-3.21-6.48⋯-2.84 0-1.030.5520.784-1.04-0.2541.31-0.526-1.141.350.416-1.830.7431.64-1.8-0.8111.79-0.3810.9782.19⋯-1.95

Real values - ⃗Re = row ( A_f; 1 ) = [1.06 -0.279 -0.942 0.809 0.555 -1.21 0.112 1.29 -0.952 -0.851 1.65 -0.148 -1.77 1.46 0.938 -2.75 1.03 4.09 -3.21 -6.48 ... -2.84]

Imaginary values - ⃗Im = row ( A_f; 2 ) = [0 -1.03 0.552 0.784 -1.04 -0.254 1.31 -0.526 -1.14 1.35 0.416 -1.83 0.743 1.64 -1.8 -0.811 1.79 -0.381 0.978 2.19 ... -1.95]

Squared magnitudes - ⃗S = ⃗Re^⊙ 2 + ⃗Im^⊙ 2 = [1.13 1.14 1.19 1.27 1.38 1.53 1.72 1.94 2.22 2.54 2.89 3.37 3.7 4.79 4.13 8.24 4.27 16.9 11.27 46.86 ... 11.87]

Fourier amplitude spectrum (FAS) - ⃗FAS = √ ⃗Sn = √ ⃗S22500 = [4.72×10^-5 4.75×10^-5 4.85×10^-5 5.01×10^-5 5.23×10^-5 5.5×10^-5 5.82×10^-5 6.2×10^-5 6.62×10^-5 7.08×10^-5 7.56×10^-5 8.16×10^-5 8.55×10^-5 9.73×10^-5 9.04×10^-5 0.000128 9.18×10^-5 0.000183 0.000149 0.000304 ... 0.000153]

Power spectral density - ⃗PSD = ⃗SFs · n = ⃗S100 Hz · 22500 = [5×10^-7 5.08×10^-7 5.3×10^-7 5.64×10^-7 6.15×10^-7 6.82×10^-7 7.63×10^-7 8.64×10^-7 9.85×10^-7 1.13×10^-6 1.28×10^-6 1.5×10^-6 1.64×10^-6 2.13×10^-6 1.84×10^-6 3.66×10^-6 1.9×10^-6 7.51×10^-6 5.01×10^-6 2.08×10^-5 ... 5.28×10^-6] s

Magnitude spectral density - ⃗ASD = √ ⃗PSD = [0.000707 0.000713 0.000728 0.000751 0.000784 0.000826 0.000874 0.000929 0.000992 0.00106 0.00113 0.00122 0.00128 0.00146 0.00136 0.00191 0.00138 0.00274 0.00224 0.00456 ... 0.0023] s^0.5

Input accelleration plot [cm/s²] vs time [s]

Fourier amplitude spectrum (FAS) [cm/s] vs frequency [Hz] in the range of 0-10 Hz

Spotted an error? Edit these examples.