sciforge_lib/maths/signal/
convolution.rs1pub fn convolve(a: &[f64], b: &[f64]) -> Vec<f64> {
2 let n = a.len() + b.len() - 1;
3 let mut out = vec![0.0; n];
4 for i in 0..a.len() {
5 for j in 0..b.len() {
6 out[i + j] += a[i] * b[j];
7 }
8 }
9 out
10}
11
12pub fn cross_correlate(a: &[f64], b: &[f64]) -> Vec<f64> {
13 let n = a.len() + b.len() - 1;
14 let mut out = vec![0.0; n];
15 for i in 0..a.len() {
16 for j in 0..b.len() {
17 out[i + b.len() - 1 - j] += a[i] * b[j];
18 }
19 }
20 out
21}
22
23pub fn autocorrelation(signal: &[f64]) -> Vec<f64> {
24 let n = signal.len();
25 let mut out = vec![0.0; n];
26 for lag in 0..n {
27 let mut s = 0.0;
28 for i in 0..n - lag {
29 s += signal[i] * signal[i + lag];
30 }
31 out[lag] = s;
32 }
33 out
34}
35
36pub fn deconvolve(signal: &[f64], kernel: &[f64]) -> Vec<f64> {
37 if kernel.is_empty() || kernel[0].abs() < 1e-30 {
38 return vec![];
39 }
40 let out_len = signal.len() - kernel.len() + 1;
41 if out_len == 0 {
42 return vec![];
43 }
44 let mut out = vec![0.0; out_len];
45 let mut remainder = signal.to_vec();
46 for i in 0..out_len {
47 out[i] = remainder[i] / kernel[0];
48 for j in 0..kernel.len() {
49 remainder[i + j] -= out[i] * kernel[j];
50 }
51 }
52 out
53}
54
55pub fn matched_filter(signal: &[f64], template: &[f64]) -> Vec<f64> {
56 let reversed: Vec<f64> = template.iter().rev().cloned().collect();
57 convolve(signal, &reversed)
58}
59
60pub fn downsample(signal: &[f64], factor: usize) -> Vec<f64> {
61 signal.iter().step_by(factor).cloned().collect()
62}
63
64pub fn upsample(signal: &[f64], factor: usize) -> Vec<f64> {
65 let mut out = vec![0.0; (signal.len() - 1) * factor + 1];
66 for (i, &s) in signal.iter().enumerate() {
67 out[i * factor] = s;
68 }
69 out
70}
71
72pub fn zero_pad(signal: &[f64], total_length: usize) -> Vec<f64> {
73 let mut out = signal.to_vec();
74 out.resize(total_length, 0.0);
75 out
76}
77
78pub fn convolve_same(a: &[f64], b: &[f64]) -> Vec<f64> {
79 let full = convolve(a, b);
80 let offset = b.len() / 2;
81 full[offset..offset + a.len()].to_vec()
82}
83
84pub fn convolve_valid(a: &[f64], b: &[f64]) -> Vec<f64> {
85 if a.len() < b.len() {
86 return vec![];
87 }
88 let out_len = a.len() - b.len() + 1;
89 let mut out = vec![0.0; out_len];
90 for i in 0..out_len {
91 let mut s = 0.0;
92 for j in 0..b.len() {
93 s += a[i + j] * b[b.len() - 1 - j];
94 }
95 out[i] = s;
96 }
97 out
98}
99
100pub fn circular_convolution(a: &[f64], b: &[f64]) -> Vec<f64> {
101 let n = a.len().max(b.len());
102 let mut ap = a.to_vec();
103 let mut bp = b.to_vec();
104 ap.resize(n, 0.0);
105 bp.resize(n, 0.0);
106 let mut out = vec![0.0; n];
107 for i in 0..n {
108 for j in 0..n {
109 out[i] += ap[j] * bp[(i + n - j) % n];
110 }
111 }
112 out
113}
114
115pub fn overlap_add(signal: &[f64], kernel: &[f64], block_size: usize) -> Vec<f64> {
116 let m = kernel.len();
117 let out_len = signal.len() + m - 1;
118 let mut output = vec![0.0; out_len];
119 let mut pos = 0;
120 while pos < signal.len() {
121 let end = (pos + block_size).min(signal.len());
122 let block = &signal[pos..end];
123 let segment = convolve(block, kernel);
124 for (i, &v) in segment.iter().enumerate() {
125 if pos + i < out_len {
126 output[pos + i] += v;
127 }
128 }
129 pos += block_size;
130 }
131 output
132}
133
134pub fn normalized_cross_correlation(a: &[f64], b: &[f64]) -> Vec<f64> {
135 let corr = cross_correlate(a, b);
136 let norm_a: f64 = a.iter().map(|x| x * x).sum::<f64>().sqrt();
137 let norm_b: f64 = b.iter().map(|x| x * x).sum::<f64>().sqrt();
138 let denom = norm_a * norm_b;
139 if denom < 1e-30 {
140 return corr;
141 }
142 corr.iter().map(|&c| c / denom).collect()
143}
144
145pub fn convolution_2d(image: &[Vec<f64>], kernel: &[Vec<f64>]) -> Vec<Vec<f64>> {
146 let ny = image.len();
147 let nx = image[0].len();
148 let ky = kernel.len();
149 let kx = kernel[0].len();
150 let hy = ky / 2;
151 let hx = kx / 2;
152 let mut out = vec![vec![0.0; nx]; ny];
153 for (j, out_j) in out.iter_mut().enumerate() {
154 for (i, out_ji) in out_j.iter_mut().enumerate() {
155 let mut s = 0.0;
156 for (kj, kernel_row) in kernel.iter().enumerate() {
157 for (ki, &kval) in kernel_row.iter().enumerate() {
158 let sj = j as i64 + kj as i64 - hy as i64;
159 let si = i as i64 + ki as i64 - hx as i64;
160 if sj >= 0 && sj < ny as i64 && si >= 0 && si < nx as i64 {
161 s += image[sj as usize][si as usize] * kval;
162 }
163 }
164 }
165 *out_ji = s;
166 }
167 }
168 out
169}
170
171pub fn wiener_deconvolution(signal: &[f64], kernel: &[f64], noise_power: f64) -> Vec<f64> {
172 let n = signal.len();
173 let (sr, si) = dft_internal(signal);
174 let mut kp = kernel.to_vec();
175 kp.resize(n, 0.0);
176 let (kr, ki) = dft_internal(&kp);
177 let mut out_re = vec![0.0; n];
178 let mut out_im = vec![0.0; n];
179 for i in 0..n {
180 let k_mag2 = kr[i] * kr[i] + ki[i] * ki[i];
181 let wiener = k_mag2 / (k_mag2 + noise_power);
182 let inv_k_re = kr[i] / (k_mag2 + 1e-30);
183 let inv_k_im = -ki[i] / (k_mag2 + 1e-30);
184 out_re[i] = wiener * (sr[i] * inv_k_re - si[i] * inv_k_im);
185 out_im[i] = wiener * (sr[i] * inv_k_im + si[i] * inv_k_re);
186 }
187 idft_internal(&out_re, &out_im)
188}
189
190fn dft_internal(signal: &[f64]) -> (Vec<f64>, Vec<f64>) {
191 let n = signal.len();
192 let mut re = vec![0.0; n];
193 let mut im = vec![0.0; n];
194 for k in 0..n {
195 for (j, &s) in signal.iter().enumerate() {
196 let angle = -2.0 * std::f64::consts::PI * k as f64 * j as f64 / n as f64;
197 re[k] += s * angle.cos();
198 im[k] += s * angle.sin();
199 }
200 }
201 (re, im)
202}
203
204fn idft_internal(re: &[f64], im: &[f64]) -> Vec<f64> {
205 let n = re.len();
206 let mut out = vec![0.0; n];
207 for (j, oj) in out.iter_mut().enumerate() {
208 for (k, (&rk, &ik)) in re.iter().zip(im.iter()).enumerate() {
209 let angle = 2.0 * std::f64::consts::PI * k as f64 * j as f64 / n as f64;
210 *oj += rk * angle.cos() - ik * angle.sin();
211 }
212 *oj /= n as f64;
213 }
214 out
215}