pub fn wavelet_transform_f32(data: &mut [f32], h: &[f32]) -> StatusExpand description
Performs a full multi-level fast wavelet transform (Ch. 27): repeatedly applies
wavelet_step_f32 to the lower half of the array, halving the active block length each
time, stopping once the block would be smaller than the filter itself (mirroring the Haar
transform’s pyramid structure).
data.len() must be a power of 2 and >= h.len().
Examples found in repository?
examples/spectral_radar_transforms.rs (line 274)
18fn main() {
19 println!("===============================================================================");
20 println!(" embedded-dsp Spectral Analysis, Multirate & Transform Pipeline ");
21 println!("===============================================================================");
22 println!();
23
24 const FS: f32 = 48000.0; // 48 kHz IF / baseband sampling rate
25 const N_SAMPLES: usize = 512;
26
27 // -----------------------------------------------------------------------------------------
28 // 1. Multitone RF / Radar Signal Synthesis with Strong Jammer
29 // -----------------------------------------------------------------------------------------
30 println!("--- 1. Multitone RF Signal Simulation with Interfering Jammer ---");
31 let mut clean_signal = [0.0f32; N_SAMPLES];
32 let mut jammer_signal = [0.0f32; N_SAMPLES];
33 let mut received_signal = [0.0f32; N_SAMPLES];
34
35 let target_f1 = 3000.0f32; // 3 kHz target tone
36 let target_f2 = 7500.0f32; // 7.5 kHz target tone
37 let jammer_freq = 1200.0f32; // 1.2 kHz strong interfering tone
38
39 for i in 0..N_SAMPLES {
40 let t = i as f32 / FS;
41 // Clean signal: Target radar return
42 clean_signal[i] = 0.6 * (2.0 * core::f32::consts::PI * target_f1 * t).sin()
43 + 0.4 * (2.0 * core::f32::consts::PI * target_f2 * t).sin();
44 // High-power narrowband interference (Jammer)
45 jammer_signal[i] = 1.8 * (2.0 * core::f32::consts::PI * jammer_freq * t).sin();
46 // Receiver noise
47 let prng =
48 ((i as u64).wrapping_mul(1664525).wrapping_add(1013904223) % 10000) as f32 / 10000.0;
49 let noise = 0.1 * (prng - 0.5);
50
51 received_signal[i] = clean_signal[i] + jammer_signal[i] + noise;
52 }
53
54 let mut raw_rms = 0.0f32;
55 rms_f32(&received_signal, &mut raw_rms);
56 println!(
57 " Generated {} samples. Total Received RMS: {:.4} (Jammer-dominated)",
58 N_SAMPLES, raw_rms
59 );
60
61 // -----------------------------------------------------------------------------------------
62 // 2. Multirate DSP: Cascaded Integrator-Comb (CIC) & Resampling
63 // -----------------------------------------------------------------------------------------
64 println!("\n--- 2. Multirate DSP: CIC Decimation / Interpolation & Linear Resampling ---");
65 // CIC Decimator: 3 stages, Decimation factor R = 4 (e.g. 48 kHz -> 12 kHz)
66 let mut cic_decimator = CicDecimator::<3>::new(4);
67 let mut cic_interpolator = CicInterpolator::<3>::new(4);
68
69 let mut decimated_stream = [0i32; 128];
70 let mut dec_count = 0;
71
72 for &sample in &received_signal {
73 let sample_i32 = (sample * 1000.0) as i32;
74 if let Some(dec_val) = cic_decimator.process_sample(sample_i32) {
75 if dec_count < decimated_stream.len() {
76 decimated_stream[dec_count] = dec_val;
77 dec_count += 1;
78 }
79 }
80 }
81 println!(
82 " CIC Decimator (R=4, 3 stages): Decimated 512 samples -> {} samples",
83 dec_count
84 );
85
86 // CIC Interpolator: Upsample by 4 back to original rate
87 let mut upsample_chunk = [0i32; 4];
88 cic_interpolator.process_sample(decimated_stream[0], &mut upsample_chunk);
89 println!(
90 " CIC Interpolator: 1 input sample -> 4 upsampled samples: {:?}",
91 upsample_chunk
92 );
93
94 // Fractional Linear Resampling (e.g., 48 kHz to 32 kHz -> ratio = 32/48 = 0.6667)
95 let mut resampled_out = [0.0f32; 341];
96 resample_linear_f32(&received_signal, &mut resampled_out, 32000.0 / 48000.0);
97 println!(
98 " Fractional Resampler (48kHz -> 32kHz): Resampled {} -> {} samples",
99 received_signal.len(),
100 resampled_out.len()
101 );
102
103 // -----------------------------------------------------------------------------------------
104 // 3. Adaptive Noise Cancellation (LMS & NLMS Filters)
105 // -----------------------------------------------------------------------------------------
106 println!("\n--- 3. Adaptive Noise Cancellation (LMS & Normalized LMS) ---");
107 // Use jammer reference signal to cancel interference from received signal
108 const LMS_TAPS: usize = 32;
109 let mut lms_coeffs = [0.0f32; LMS_TAPS];
110 let mut lms_state = [0.0f32; LMS_TAPS];
111 let mut lms_filter =
112 LmsInstanceF32::init(LMS_TAPS as u16, &mut lms_coeffs, &mut lms_state, 0.005);
113
114 let mut lms_cancelled_out = [0.0f32; N_SAMPLES];
115 let mut lms_error = [0.0f32; N_SAMPLES];
116
117 // Reference signal x = jammer reference, Desired d = received signal (target + jammer + noise)
118 // Error output e = d - y = target + noise (jammer cancelled!)
119 lms_f32(
120 &mut lms_filter,
121 &jammer_signal,
122 &received_signal,
123 &mut lms_cancelled_out,
124 &mut lms_error,
125 );
126
127 let mut err_rms_initial = 0.0f32;
128 let mut err_rms_converged = 0.0f32;
129 rms_f32(&lms_error[..64], &mut err_rms_initial);
130 rms_f32(&lms_error[N_SAMPLES - 64..], &mut err_rms_converged);
131
132 println!(" LMS Adaptive Filter (32 taps, mu=0.005):");
133 println!(" • Initial Error RMS : {:.4}", err_rms_initial);
134 println!(
135 " • Converged Error RMS : {:.4} (Interference cancelled, target recovered!)",
136 err_rms_converged
137 );
138
139 // -----------------------------------------------------------------------------------------
140 // 4. Windowing Comparison (Spectral Leakage Suppression)
141 // -----------------------------------------------------------------------------------------
142 println!("\n--- 4. Windowing Comparison: Mainlobe & Sidelobe Properties ---");
143 const WIN_LEN: usize = 64;
144 let win_rect = [1.0f32; WIN_LEN];
145 let mut win_hanning = [0.0f32; WIN_LEN];
146 let mut win_hamming = [0.0f32; WIN_LEN];
147 let mut win_blackman_harris = [0.0f32; WIN_LEN];
148 let mut win_flat_top = [0.0f32; WIN_LEN];
149
150 hanning_f32(&mut win_hanning);
151 hamming_f32(&mut win_hamming);
152 blackman_harris_f32(&mut win_blackman_harris);
153 flattop_f32(&mut win_flat_top);
154
155 println!(" Sample Window Values at Midpoint (n=32):");
156 println!(" • Rectangular : {:.4}", win_rect[32]);
157 println!(" • Hanning : {:.4}", win_hanning[32]);
158 println!(" • Hamming : {:.4}", win_hamming[32]);
159 println!(" • Blackman-Harris: {:.4}", win_blackman_harris[32]);
160 println!(" • Flat-Top : {:.4}", win_flat_top[32]);
161
162 // -----------------------------------------------------------------------------------------
163 // 5. High-Resolution Spectral Analysis (FFT & Welch PSD)
164 // -----------------------------------------------------------------------------------------
165 println!("\n--- 5. Spectral Analysis: In-Place Complex FFT & Welch's PSD ---");
166 // 256-point Complex FFT on Converged LMS Cleaned Signal
167 const FFT_SIZE: usize = 256;
168 let mut fft_buffer = [0.0f32; 2 * FFT_SIZE];
169 for i in 0..FFT_SIZE {
170 fft_buffer[2 * i] = lms_error[N_SAMPLES - FFT_SIZE + i] * win_hamming[i % WIN_LEN];
171 fft_buffer[2 * i + 1] = 0.0;
172 }
173
174 cfft_f32(&mut fft_buffer, FFT_SIZE, 0, 1);
175
176 // Find top 2 spectral peaks
177 let mut peak1_mag = 0.0f32;
178 let mut peak1_bin = 0usize;
179 let mut peak2_mag = 0.0f32;
180 let mut peak2_bin = 0usize;
181
182 for k in 1..FFT_SIZE / 2 {
183 let re = fft_buffer[2 * k];
184 let im = fft_buffer[2 * k + 1];
185 let mag = (re * re + im * im).sqrt();
186 if mag > peak1_mag {
187 peak2_mag = peak1_mag;
188 peak2_bin = peak1_bin;
189 peak1_mag = mag;
190 peak1_bin = k;
191 } else if mag > peak2_mag {
192 peak2_mag = mag;
193 peak2_bin = k;
194 }
195 }
196
197 let bin_to_hz = FS / FFT_SIZE as f32;
198 println!(" 256-Point FFT Detected Spectral Peaks:");
199 println!(
200 " • Peak 1: Bin {:>3} ({:>6.1} Hz) -> Magnitude: {:>7.2}",
201 peak1_bin,
202 peak1_bin as f32 * bin_to_hz,
203 peak1_mag
204 );
205 println!(
206 " • Peak 2: Bin {:>3} ({:>6.1} Hz) -> Magnitude: {:>7.2}",
207 peak2_bin,
208 peak2_bin as f32 * bin_to_hz,
209 peak2_mag
210 );
211
212 // Welch's Method Power Spectral Density (PSD)
213 let mut psd_out = [0.0f32; 64];
214 welch_psd_f32(
215 &received_signal,
216 &mut psd_out,
217 128, // Segment length
218 64, // Overlap
219 FS, // Sampling rate
220 WelchWindow::BlackmanHarris,
221 true, // Return in dB/Hz
222 );
223 println!(" Welch's PSD (Averaged Periodogram in dB/Hz):");
224 println!(" • PSD[Bin 4] (375 Hz) : {:>6.1} dB/Hz", psd_out[4]);
225 println!(
226 " • PSD[Bin 13] (1200 Hz): {:>6.1} dB/Hz (Jammer peak)",
227 psd_out[13]
228 );
229 println!(
230 " • PSD[Bin 32] (3000 Hz): {:>6.1} dB/Hz (Target 1 peak)",
231 psd_out[32]
232 );
233
234 // -----------------------------------------------------------------------------------------
235 // 6. Advanced DSP Transforms (FWHT, DCT-IV, Hartley, Wavelets)
236 // -----------------------------------------------------------------------------------------
237 println!("\n--- 6. Advanced DSP Transforms: FWHT, DCT-IV, Hartley & Wavelets ---");
238
239 // A. Fast Walsh-Hadamard Transform (FWHT) for CDMA / Walsh codes
240 let mut walsh_data = [1.0f32, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0];
241 let fwht_status = fwht_f32(&mut walsh_data);
242 println!(
243 " Fast Walsh-Hadamard Transform (8-point): Status={:?}, Out={:?}",
244 fwht_status, walsh_data
245 );
246 ifwht_f32(&mut walsh_data);
247 println!(
248 " Inverse FWHT (reconstructed exact signal): {:?}",
249 walsh_data
250 );
251
252 // B. Discrete Cosine Transform IV (DCT-IV)
253 let dct_in = [1.0f32, 2.0, 3.0, 4.0];
254 let mut dct_out = [0.0f32; 4];
255 dct4_f32(&dct_in, &mut dct_out, 4);
256 println!(" DCT-IV (4-point): In={:?} -> Out={:?}", dct_in, dct_out);
257
258 // C. Hartley Transform (Real-in, Real-out, Self-Inverse Transform)
259 let mut hartley_buf = [1.0f32, 2.0, 3.0, 4.0];
260 let h_status = hartley_transform_f32(&mut hartley_buf);
261 println!(
262 " Hartley Transform: Status={:?}, F_H={:?}",
263 h_status, hartley_buf
264 );
265 // Applying Hartley twice recovers the original scaled signal
266 hartley_transform_f32(&mut hartley_buf);
267 println!(" Hartley Roundtrip (Self-Inverse): {:?}", hartley_buf);
268
269 // D. Multi-Resolution Daubechies-4 Discrete Wavelet Transform (DWT)
270 let mut wavelet_signal = [0.0f32; 16];
271 wavelet_signal[7] = 10.0; // High-frequency transient spike in middle of frame
272 let orig_wavelet = wavelet_signal;
273
274 let dwt_status = wavelet_transform_f32(&mut wavelet_signal, &DAUBECHIES_4);
275 println!(" Daubechies-4 DWT (16-point Pyramid Decomposition):");
276 println!(" • DWT Status : {:?}", dwt_status);
277 println!(
278 " • Detail/Scaling Coefficients: {:?}",
279 &wavelet_signal[..8]
280 );
281
282 let idwt_status = inverse_wavelet_transform_f32(&mut wavelet_signal, &DAUBECHIES_4);
283 let mut diff = 0.0f32;
284 for i in 0..16 {
285 diff += (wavelet_signal[i] - orig_wavelet[i]).abs();
286 }
287 println!(
288 " • Inverse DWT Status: {:?}, Reconstruction Error: {:.2e}",
289 idwt_status, diff
290 );
291
292 println!();
293 println!("===============================================================================");
294 println!(" Spectral & Transforms Pipeline Execution Complete! ");
295 println!("===============================================================================");
296}