Waves In Ice Observing System (WIIOS) Development
During research into wave-ice interactions, I discovered a significant disconnect: while theoretical models were becoming increasingly sophisticated, we lacked the observational data needed to validate them. This gap, combined with predictions of intensifying wave activity in the Southern Ocean, motivated me to transition from theoretical modeling to hands-on data collection.
Developing Innovative Technology
The harsh polar environment presented unique challenges to obtaining observational data: massive waves, crushing ice, hurricane-force winds, and sub-zero temperatures demanded exceptional equipment durability. In collaboration with P.A.S. software and electrical engineers, we developed custom wave-in-ice sensors. This project benefited from the expertise of Martin Doble, whose insights were crucial for hull design and wave data analysis. While P.A.S. focused on hardware, software, and GUI development, I specialized in implementing onboard wave analysis algorithms to ensure precise data collection. See Kohout et al. (2015) and Kohout et al. (2020) for a detailed discussion of the system. Below, I outline the core components of the wave motion processing, accompanied by relevant code snippets and GitHub links.
Real-Time Spectral Analysis for Wave Motion Processing in Ice
As part of the Waves in Ice Observation System (WIIOS) project, I developed a C-based signal processing system to analyze wave motion data in sea ice environments. This project was a unique challenge, as my formal programming training was limited to Fortran basics during my mathematics degree. I independently learned C to build this system from the ground up. Due to the constrained memory on the sensors, only essential libraries were included, requiring me to write all algorithms from scratch. As a result, the code deviates from conventional formatting standards. Please note that this software was initially developed for personal use and not intended for public sharing. I am currently working on refining the code with AI assistance to make it accessible in a polished, shareable format.
Data Acquisition & Preprocessing
[GitHub Links: main.c, sc_misc.c, signal_conditioning.c]
- Captured 640-second burst samples at 64 Hz using IMU and Kistler accelerometers.
- Designed a digital filtering pipeline, including:
- A low-pass 2nd-order Butterworth filter (0.5 Hz cutoff).
- Downsampling to 2 Hz for efficient processing.
Code Snippet: Downsampling Implementation
/******************************************************************************
* Function: decimate
* Purpose: Decimates (downsamples) a signal after applying an anti-aliasing
* Butterworth filter
*
* Parameters:
* in[] - Input signal array
* out[] - Output decimated signal array
* sr1 - Input sampling rate (Hz)
* sr2 - Output sampling rate (Hz)
* n1 - Length of input array
* n2 - Length of output array
*
* Returns:
* void
*
* Notes:
* - Implements 2nd order Butterworth filter
* - Supports two specific sampling rate combinations:
* 1. 8 Hz input (sr1=0.125) with 0.5 Hz cutoff
* 2. 64 Hz input (sr1=0.015625) with 2 Hz cutoff
* - Requires integer ratio between input and output sampling rates
* - Includes bounds checking and error logging
*
* Filter Coefficients:
* 8 Hz sampling: a = [1.0, -1.4542, 0.5741]
* b = [0.0300, 0.0599, 0.0300]
* 64 Hz sampling: a = [1.0, -1.7238, 0.7575]
* b = [0.0084, 0.0169, 0.0084]
*
*****************************************************************************/
void decimate(long double *in, long double *out, float sr1, float sr2, int n1, int n2) {
logger("Start decimate",0.0);
int i, j, ii, n, jj;
int step;
long double out1[n1];
long double a[3],b[3];
logger(" n1",n1);
// intialize out
for (i=0;i<n2;i++) out[i] = 0.0;
// butter worth filter with 8 Hz sample rate and 0.5 Hz cut off
if (sr1==0.125) {
logger(" Butter worth filter with 8 Hz sample rate, 0.5 Hz cut off",0.0);
a[0] = 1.000000000000000;
a[1] = -1.454243586251585;
a[2] = 0.574061915083955;
b[0] = 0.029954582208092;
b[1] = 0.059909164416185;
b[2] = 0.029954582208092;
}
// butter worth filter with 64 Hz sample rate and 2 Hz cut off
if (sr1==0.015625) {
logger(" Butter worth filter with 64 Hz sample rate, 2 Hz cut off",0.0);
a[0] = 1.000000000000000;
a[1] = -1.723776172762509;
a[2] = 0.757546944478829;
b[0] = 0.008442692929080;
b[1] = 0.016885385858160;
b[2] = 0.008442692929080;
}
if (!floorf(sr2/sr1==sr2/sr1)) {
debug("Error: step a non integer",0.0);
return;
}
out1[0] = in[0];
out1[1] = in[1];
for (i=2; i<n1; i++){
out1[i] = b[0]*in[i] + b[1]*in[i-1] + b[2]*in[i-2]
- a[1]*out1[i-1] - a[2]*out1[i-2];
}
step = sr2/sr1;
logger(" step",step);
for (i=0, j=0; i<n1; i+=step,j++){
if (j>=n2) {
fprintf(stderr,"ERROR: out of bounds in decimate. Exiting.\n");
exit(0);
}
out[j] = out1[i];
}
logger("End decimate",0.0);
return;
}
Spectral Analysis & Directional Estimation
[GitHub Links: process.c, nr_four.c, direction.c]
Implemented Welch's method for estimating the total energy distribution across frequencies (power spectral density - PSD).
- 256-second segments with 50% overlap
- 10% cosine windowing
- Detrending algorithms for each segment
- Spectral moment calculations (derived from PSD)
Using cross-spectral analysis between vertical and horizontal motion data, the system estimates wave direction and spread. Multiple methods were explored and implemented, including approaches similar to the First-order Extended Maximum Likelihood Method (FEM) and the Maximum Entropy Method (MEM) (or Kuik et al.'s method).
- Calculated cross-spectral densities between vertical acceleration and horizontal motions (pitch/roll).
- Derived directional parameters per frequency band, including:
- Peak Wave Direction (direction of maximum energy)
- Mean Wave Direction
- Directional Spread (measure of how focused or broad the wave directions are)
- Directional quality parameters (e.g., ratio of directional components)
- Integrated results over relevant frequency ranges to obtain overall directional statistics.
Code Snippet: Welch's Method PSD Estimation with Cosine Tapering
/******************************************************************************
* Function: spctrm_psd
* Purpose: Calculates Power Spectral Density using Welch's method with
* overlapping segments and 10% cosine tapering
*
* Parameters:
* disp[] - Input displacement time series data
* p[] - Output array for power spectral density values
* m - Half the segment length (full segment = 2m)
* k - Number of segments to average (determines overlap)
*
* Returns:
* 0 on success
*
* Notes:
* - Applies detrending and demeaning to each segment
* - Uses 10% cosine taper window
* - Includes overlap processing for improved spectral estimates
* - Normalizes output by window power and number of segments
*
*****************************************************************************/
int spctrm_psd(long double disp[], long double p[], int m, int k) {
FILE *id;
int i=0,ii=0,j=0,count=0;
int mm=2*m;
int gamma;
int cl = 1;
float r=0.0;
long double win[mm], sumw=0;
long double dp[mm], psd[2*k][m+1];
long double alpha, beta, t1;
long double w1[mm];
// Intialize
for (i=0;i<mm;i++) dp[i] = 0.0;
for (i=0;i<m+1;i++) {
for(j=0;j<2*k;j++) psd[j][i] = 0.0;
p[i] = 0.0;
}
// define 10% cosine taper
gamma = 5; t1 = m-m/gamma; alpha = 2.0*3.14159265*gamma/mm;
if (gamma % 2 == 0) beta = 3.14159265; else beta = 0;
for (i=0;i<m;i++) {
if (i<t1) win[i] = 1;
else win[i] = 0.5*(1 + cos(alpha*i + beta));
}
for (i=m;i<mm;i++) win[i] = win[i-m];
for (i=0;i<m;i++) win[i] = win[mm-1-i];
for (i=0;i<mm;i++) sumw += pow(win[i],2.0);
for (i=0;i<mm;i++) if(isnan(win[i])) debug("nan",i);
for (j=0,count=0; j<2*k; j++, count+=m) {
for (i=0;i<mm;i++) {
dp[i] = disp[count+i];
demean(dp,i);
trend_removal(dp,mm);
demean(dp,i);
w1[i] = win[i]*dp[i];
}
realft(w1,mm,1);
psd[j][0] = pow(w1[0],2); psd[j][m] = pow(w1[1],2);
for (i=2,ii=1; i<mm-1; i+=2,ii++) {
if (ii>=m+1) {
fprintf(stderr,"ERROR: out of bounds in spctrm_psd. Exiting.\n");
exit(0);
}
if (i>=2*m-1) {
fprintf(stderr,"ERROR: out of bounds in spctrm_psd. Exiting.\n");
exit(0);
}
psd[j][ii] = 2*(pow(w1[i],2) + pow(w1[i+1],2));
}
}
for (i=0; i<m+1; i++) {
p[i]=0.0;
for (j=0; j<2*k; j++) p[i]+=psd[j][i];
p[i] /= (2.0*k*mm*sumw);
}
return(0);
}
Code Snippet: Wave Directional Analysis Using First-Order Extended Maximum Likelihood Method
/******************************************************************************
* Function: wave_direction_FEM_2hz_accel
* Purpose: Estimates wave direction using FEM (First-order Extended Maximum
* likelihood Method) from 2Hz accelerometer and motion sensor data
*
* Parameters:
* accel[] - Vertical acceleration time series
* pitch[] - Pitch motion time series
* roll[] - Roll motion time series
* return_data - Struct containing output parameters including:
* - direction_full[] (Full directional spectrum)
* - peak_direction (Direction at spectral peak)
* - direction (Mean direction)
* - spread (Directional spread)
* - ratio (Directional quality parameter)
* - hs_dir (Significant wave height)
*
* Returns:
* 0 on success
*
* Notes:
* - Uses cross-spectral analysis between vertical and horizontal motions
* - Frequency range defined by rspns_f7 to rspns_f8
* - Outputs directional parameters in degrees (0° = waves toward East)
* - Implements circular statistics for directional averaging
* - Includes spectral analysis using Welch's method
*
* References:
* - Extended Maximum Likelihood Method for wave direction analysis
* - Circular statistics for directional spread calculation
*
* Output Files:
* - roll_spectra.out - Roll motion spectra
* - pitch_spectra.out - Pitch motion spectra
* - accel_spectra.out - Acceleration spectra
* - direction_full.out - Full directional spectrum
* - spread.out - Directional spread
*
*****************************************************************************/
int wave_direction_FEM_2hz_accel(long double *accel, long double *pitch, long double *roll, struct RETDATA *return_data){
logger("Start wave_direction_FEM_2hz_accel",0);
int i,j,peak_id;
int w=m, ww=2*w, kk=2;
int st, ed;
int send_dir[w+1];
long double df = 1/(ww*dt);
long double c12[w+1], q12[w+1];
long double c13[w+1], q13[w+1];
long double c11[w+1], c22[w+1], c33[w+1];
long double kd[w+1];
long double fft_accel[ww], fft_ns[ww], fft_ew[ww];
long double a1[w+1], b1[w+1];
long double a2[w+1], b2[w+1];
long double a3[w+1], b3[w+1];
long double alpha12[w+1], alpha13[w+1];
long double aa1[w+1], bb1[w+1], dir[w+1], sprd[w+1], ui[w+1];
long double aa1s[w+1], bb1s[w+1];
long double peak = 0, omega, mean_dir;
long double Hs, c11_sum, dir_mean, ui_mean=0;
long double ff[w+1];
// Define ff
for (i=0;i<w+1;i++) ff[i] = (i/(2.0*m*dt));
// Initialize variables
for (i=0;i<w+1;i++){
c11[i] = 0.0;
c22[i] = 0.0;
c33[i] = 0.0;
c12[i] = 0.0;
q12[i] = 0.0;
c13[i] = 0.0;
q13[i] = 0.0;
a1[i] = 0.0;
a2[i] = 0.0;
a3[i] = 0.0;
b1[i] = 0.0;
b2[i] = 0.0;
b3[i] = 0.0;
alpha12[i] = 0.0;
alpha13[i] = 0.0;
aa1[i] = 0.0;
bb1[i] = 0.0;
dir[i] = 0.0;
sprd[i] = 0.0;
send_dir[i] = 0;
}
for (i=0;i<ww;i++){
fft_accel[i] = 0.0;
fft_ns[i] = 0.0;
fft_ew[i] = 0.0;
}
spctrm(roll,fft_ew,w,kk,n_spec,(float)dt);
spctrm(pitch,fft_ns,w,kk,n_spec,(float)dt);
spctrm(accel,fft_accel,w,kk,n_spec,(float)dt);
print2out(fft_ew,ww,"out/roll_spectra.out",ido_rs);
print2out(fft_ns,ww,"out/pitch_spectra.out",ido_ps);
print2out(fft_accel,ww,"out/accel_spectra.out",ido_hs);
// auto / cross correlations
for (i=0,j=0; i<ww; i+=2,j++){
if (j>=w+1) {
fprintf(stderr,"ERROR: out of bounds in wave_direction_FEM_2hz. Exiting.\n");
exit(0);
}
a1[j] = fft_accel[i];
a2[j] = fft_ns[i];
a3[j] = fft_ew[i];
}
a1[w] = fft_accel[ww-1];
a2[w] = fft_ns[ww-1];
a3[w] = fft_ew[ww-1];
b1[0] = 0;
b1[1] = 0;
b2[0] = 0;
b2[1] = 0;
b3[0] = 0;
b3[1] = 0;
for (i=3,j=2; i<ww; i+=2,j++){
if (j>=w+1) {
fprintf(stderr,"ERROR: out of bounds in wave_direction_FEM_2hz. Exiting.\n");
exit(0);
}
b1[j] = fft_accel[i];
b2[j] = fft_ns[i];
b3[j] = fft_ew[i];
}
b1[w] = 0;
b2[w] = 0;
b3[w] = 0;
// Find range for directional analysis
for (i = 0; i<w+1; i++) {
if (ff[i] > rspns_f7) {
st = i;
break;
}
}
for (i = 0; i<w+1; i++) {
if (ff[i] > rspns_f8) {
ed = i;
break;
}
}
c11_sum = 0;
dir_mean = 0;
Hs = 0;
for (i = 1; i<m+1; i++){
omega = (2.0*pi*ff[i])*(2.0*pi*ff[i]);
c12[i] = a1[i]*a2[i]/omega + b1[i]*b2[i]/omega;
q12[i] = a1[i]*b2[i]/omega - b1[i]*a2[i]/omega;
c13[i] = a1[i]*a3[i]/omega + b1[i]*b3[i]/omega;
q13[i] = a1[i]*b3[i]/omega - b1[i]*a3[i]/omega;
c11[i] = a1[i]*a1[i]/(omega*omega) + b1[i]*b1[i]/(omega*omega);
c22[i] = (a2[i]*a2[i] + b2[i]*b2[i]);
c33[i] = (a3[i]*a3[i] + b3[i]*b3[i]);
alpha12[i] = atan2(c12[i],q12[i]);
alpha13[i] = atan2(c13[i],q13[i]);
aa1[i] = (c12[i]*sin(alpha12[i]) + q12[i]*cos(alpha12[i]))/sqrt(c11[i]*(c22[i]+c33[i]));
bb1[i] = (c13[i]*sin(alpha13[i]) + q13[i]*cos(alpha13[i]))/sqrt(c11[i]*(c22[i]+c33[i]));
kd[i] = sqrt((c22[i]+c33[i])/c11[i]);
aa1s[i] = q12[i]/(kd[i]*c11[i]);
bb1s[i] = q13[i]/(kd[i]*c11[i]);
// note that 0 deg is for waves headed towards positive x (EAST, right hand system)
dir[i] = atan2(bb1[i],aa1[i]);
ui[i] = sqrt(aa1s[i]*aa1s[i] + bb1s[i]*bb1s[i]);
sprd[i] = sqrt(2.0*(1.0-ui[i])); //circular spread (Degrees) from Doble
return_data->direction_full[i] = round(dir[i]*1800.0/pi);
if ((i >= st) && (i < ed)) {
c11_sum += c11[i];
dir_mean += dir[i];
ui_mean += ui[i];
if (peak<=c11[i]) {
peak = c11[i];
peak_id = i;
}
}
}
return_data->peak_direction = dir[peak_id];
return_data->direction = dir_mean/(ed-st);
return_data->spread = sprd[peak_id];
return_data->ratio = ui_mean/(ed-st);
//return_data->ratio =sqrt(c11[peak_id]/(c22[peak_id]/kd[peak_id] + c33[peak_id]/kd[peak_id]));
// Significant wave height
return_data->hs_dir = 4.0*sqrt(c11_sum);
print2out(c11,w+1,"out/c11.out",ido_od);
print2out(dir,w+1,"out/direction_full.out",ido_od);
print2out(sprd,w+1,"out/spread.out",ido_od);
print2out(aa1,w+1,"out/aa1.out",ido_od);
print2out(bb1,w+1,"out/bb1.out",ido_od);
logger("End wave_direction_FEM_2hz",0);
return(0);
}
Automated Quality Control System
[GitHub Link: qc.c]
- Developed validation mechanisms, including statistical tests.
- Implemented spike detection to identify anomalies.
- Monitored consecutive data changes to ensure consistency.
Code Snippet: Accelerometer Data Quality Control and Signal Processing
/******************************************************************************
* Function: quality_control_accel
* Purpose: Performs quality control checks on accelerometer data including
* spike detection, flat signal detection, and statistical validation
*
* Parameters:
* data_raw - Struct containing raw accelerometer data:
* - accel_raw_kist (Kistler accelerometer)
* - accel_raw_x/y/z (IMU accelerometer axes)
* return_data - Struct containing quality flags and metrics:
* - qflg_kist (Kistler quality flag)
* - qflg_accel (IMU quality flag)
* - qflg_pkist (Kistler percentage flag)
*
* Returns:
* void
*
* Quality Checks:
* 1. Flat Signal Detection
* - Identifies unresponsive periods
* - Threshold: 0.5
* - Applies patch correction
*
* 2. Spike Detection
* - Kistler: 10 iterations with nstd_kist threshold
* - IMU: Single pass with nstd_imu threshold
* - Applies patch correction after detection
*
* 3. Statistical Validation
* - Performs basic statistical tests
* - Sets quality flags based on results
*
* Output Metrics:
* - Percentage of compromised data
* - Counts of detected spikes and flat periods
* - Separate quality flags for Kistler and IMU data
*
* Notes:
* - Implements different thresholds for Kistler and IMU
* - Patches detected issues in-place
* - Generates detailed quality control log if PRINT_OUTPUT defined
*
*****************************************************************************/
void quality_control_accel(struct DATARAW *data_raw, struct RETDATA *return_data){
logger("Start quality_control_accel",0.0);
int i, r;
int t1 = 0;
int n_flat_k=0,n_flat_x=0, n_flat_y=0, n_flat_z=0;
int max_flat_k=0, max_flat_x=0, max_flat_y=0, max_flat_z=0;
int n_spike_k=0, n_spike_x=0, n_spike_y=0, n_spike_z=0;
int max_spike_k=0, max_spike_x=0, max_spike_y=0, max_spike_z=0;
int qflg_k=0,qflg_x=0,qflg_y=0,qflg_z=0;
int temp_n, temp_max;
// Flat
qflg_k = flat(data_raw->accel_raw_kist,&n_flat_k,&max_flat_k,0.5,n_raw);
if (qflg_k == 1) return_data->qflg_kist = 1;
patch(data_raw->accel_raw_kist,n_raw);
qflg_x = flat(data_raw->accel_raw_x,&n_flat_x,&max_flat_x,0.5,n_raw);
patch(data_raw->accel_raw_x,n_raw);
qflg_y = flat(data_raw->accel_raw_y,&n_flat_y,&max_flat_y,0.5,n_raw);
patch(data_raw->accel_raw_y,n_raw);
qflg_z = flat(data_raw->accel_raw_z,&n_flat_z,&max_flat_z,0.5,n_raw);
patch(data_raw->accel_raw_z,n_raw);
// Despike
// run despike on Kistler several times to remove all spikes
for (i=0;i<10;i++){
despike(data_raw->accel_raw_kist,&temp_n,&temp_max,n_raw,nstd_kist);
patch(data_raw->accel_raw_kist,n_raw);
n_spike_k = n_spike_k + temp_n;
max_spike_k = max_spike_k + temp_max;
}
despike(data_raw->accel_raw_x,&n_spike_x,&max_spike_x,n_raw,nstd_imu);
patch(data_raw->accel_raw_x,n_raw);
despike(data_raw->accel_raw_y,&n_spike_y,&max_spike_y,n_raw,nstd_imu);
patch(data_raw->accel_raw_y,n_raw);
despike(data_raw->accel_raw_z,&n_spike_z,&max_spike_z,n_raw,nstd_imu);
patch(data_raw->accel_raw_z,n_raw);
// Test stats
if (return_data->qflg_kist == 0) {
return_data->qflg_kist = test_stats(data_raw->accel_raw_kist,n_raw);
}
r = test_stats(data_raw->accel_raw_x,n_raw);
if (r == 1) return_data->qflg_accel = 100;
r = test_stats(data_raw->accel_raw_y,n_raw);
if (r == 1) return_data->qflg_accel = 100;
r = test_stats(data_raw->accel_raw_z,n_raw);
if (r == 1) return_data->qflg_accel = 100;
// Return flags
if (return_data->qflg_accel <100) {
return_data->qflg_accel = (int)100*(n_flat_x + n_spike_x + n_flat_y + n_spike_y + n_flat_z + n_spike_z)/(n_raw*3);
}
return_data->qflg_pkist = (int)100*(n_flat_k + n_spike_k)/(n_raw);
// Print
#ifdef PRINT_OUTPUT
fprintf(ido_oqc,"ACCEL QUALITY FLAGS\n");
fprintf(ido_oqc,"Quality flag x = %d\n", qflg_x);
fprintf(ido_oqc,"Quality flag y = %d\n", qflg_y);
fprintf(ido_oqc,"Quality flag z = %d\n", qflg_z);
fprintf(ido_oqc,"Number of unresponsive x = %d\n",n_flat_x);
fprintf(ido_oqc,"Number of unresponsive y = %d\n",n_flat_y);
fprintf(ido_oqc,"Number of unresponsive z = %d\n",n_flat_z);
fprintf(ido_oqc,"Number of spikes x = %d\n",n_spike_x);///(long double)n_raw);
fprintf(ido_oqc,"Number of spikes y = %d\n",n_spike_y);///(long double)n_raw);
fprintf(ido_oqc,"Number of spikes z = %d\n",n_spike_z);///(long double)n_raw);
#endif
logger("End quality_control_accel",0.0);
return ;
}
Real-Time Processing Constraints
- Utilized a dual-core Edison processor for computations.
- Managed operations within a 1 GB RAM limitation.
- Optimized data storage on a 32 GB SD card.
Technical Skills Highlighted
- Advanced digital signal processing techniques.
- Real-time programming for embedded systems.
- Statistical analysis and data validation methods.
- Sensor integration and memory optimization.
This project demonstrates a blend of advanced programming proficiency and deep mathematical understanding of spectral analysis and signal processing, tailored to the unique challenges of studying wave dynamics in Antarctic sea ice.