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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
use super::{CPSSettings, Measurement, error::*, load_store_cps::*};
use crate::file_attr_names::MEASURED_INPUT_DATASET_NAME;
use crate::ps::{CPSResult, getFreq};
use crate::{
measurement::{
cpsthread::{CPSComputeStatus, CPSThread},
meas_helper::ConvertedBlockIter,
},
*,
};
use hdf5_metno::File;
use snafu::prelude::*;
type Result<T> = std::result::Result<T, MeasurementError>;
impl Measurement {
/// Compute / load Cross Power Spectra for measurement, based on settings
pub fn CPS(
&mut self,
settings: &CPSSettings,
fw: FreqWeighting,
channels: Option<&[usize]>,
) -> Result<CPSResult> {
self.joinFinishedCPSThreads(Some(settings));
// Compute / restore Cross Power Spectra
let cps = {
let file = self.open_file(false)?;
let cps = load_CPS(&file, settings)?;
drop(file);
if let Some(cps) = cps {
cps
} else {
let cps = self.compute_CPS(settings)?;
let cps = cps.with_context(|| {
let istart = settings.istart.unwrap_or(0);
let istop = settings.istop.unwrap_or(self.nframes);
NFFTTooLargeSnafu {
nfft: settings.nfft,
nsamples: istop - istart,
}
})?;
// We computed new spectra, now we cache them in the file
let file = self.open_file(true)?;
Self::storeCPSResultInCache(&file, settings, &cps)?;
cps
}
};
// Channel selection
let (mut cps, nchannels) = if let Some(channels) = channels {
ensure!(
!channels.is_empty(),
LogicSnafu {
name: self.name(),
message: "No channels provided"
}
);
let highest_channel = channels
.iter()
.copied()
.max()
.expect("No channels provided, but check passed?");
ensure!(
highest_channel < self.nchannels(),
ChannelIdxOutOfBoundsSnafu {
channel_idx: highest_channel,
max_channels: self.nchannels()
}
);
// Copy over only channels that are selected.
let mut ch_sel = Array3::zeros((cps.shape()[0], channels.len(), channels.len()));
for i in 0..channels.len() {
for j in 0..channels.len() {
ch_sel.slice_mut(s![.., i, j]).assign(&cps.slice(s![
..,
channels[i],
channels[j]
]));
}
}
let sens = self.sensitivities(Some(channels));
let mut res = CPSResult::new(ch_sel);
res.apply_sensitivities(&sens);
(res, channels.len())
} else {
let mut cps = CPSResult::new(cps);
let sens = self.sensitivities(None);
cps.apply_sensitivities(&sens);
(cps, self.nchannels())
};
// Apply frequency weighting, if not Z-weighting
if !matches!(fw, FreqWeighting::Z) {
let freq = getFreq(self.samplerate(), settings.nfft);
assert!(freq.len() == cps.inner().shape()[0]);
let freq_slice = freq.as_slice().expect("Slicing 1D should work");
let sq_weight = fw
.powerweight(freq_slice)
.into_shape_with_order((freq_slice.len(), 1, 1))
.expect("Cannot create shape");
let sq_weight = sq_weight
.broadcast((freq_slice.len(), nchannels, nchannels))
.expect("Broadcasting should work");
Zip::from(cps.inner_mut())
.and(sq_weight)
.for_each(|c, w| *c *= w);
}
Ok(cps)
}
/// Clear the cached CPS data from the measurement file. After calling
/// this, any subsequent call to [Measurement::CPS] will recompute the
/// spectra from the raw audio data.
pub fn clearCPSCache(&mut self) -> Result<()> {
let file = self.open_file(true)?;
clear_all_cached_cps(&file)?;
Ok(())
}
/// Cleanup all running CPS threads. Do not process there data.
pub fn killRunningCPSThreads(&mut self) {
self.cps_threads.clear();
}
/// Add a new CPS thread to the measurement, that is computing Cross Power
/// Spectra for the data already stored in the measurement file. It might be
/// done, in which case the data is stored in the measurement file. When CPS data is requested,
pub fn add_cpsthread(&mut self, cpsthread: CPSThread) {
self.cps_threads.push(cpsthread);
}
// -- Private methods below
/// Compute power spectra, call this in case CPS is not cached.
///
/// # Arguments - `settings`: CPSSettings struct containing settings for
/// power spectrum computation
///
/// # Returns - `Result<Option<Array3<Cflt>>>`: Result containing the power
/// spectra or None in case too little data is available
fn compute_CPS(&self, settings: &CPSSettings) -> Result<Option<Array3<Cflt>>> {
// let freq = getFreq(self.samplerate, settings.nfft);
let istart_istop = if let Some(istart) = settings.istart {
if let Some(istop) = settings.istop {
Some((istart, istop))
} else {
Some((istart, self.nframes))
}
} else {
None
};
let mut cpsthread = CPSThread::new(self.samplerate(), *settings, Some(5));
// Note: no sensitivities are applied here, to match the CPS that is
// computed during the recording (which uses uncalibrated data).
let blockiter = ConvertedBlockIter::new(
self,
MEASURED_INPUT_DATASET_NAME,
istart_istop,
None,
None,
None,
)?;
for block in blockiter {
cpsthread.push(block.into());
}
let res = cpsthread.join();
Ok(res.1)
}
/// Store Cross Power Spectra in cache file. Small helper function
fn storeCPSResultInCache(
file: &File,
settings: &CPSSettings,
cps: &Array3<Cflt>,
) -> Result<()> {
store_CPS(file, settings, cps)?;
Ok(())
}
/// Join running CPS Threads that are done, putting the result in the
/// measurement file.
///
/// # Args:
///
/// - `waitFor`: If given, wait for a CPS thread with these settings to be
/// done before returning
fn joinFinishedCPSThreads(&mut self, waitFor: Option<&CPSSettings>) {
let nthreads = self.cps_threads.len();
// Take ownership of the threads, leaving an empty vector in its place.
let threads = std::mem::replace(&mut self.cps_threads, Vec::with_capacity(nthreads));
if threads.is_empty() {
return;
}
match self.open_file(true) {
Ok(file) => {
threads.into_iter().for_each(|mut thread| {
let settings = *thread.getSettings();
match thread.get_status() {
CPSComputeStatus::Waiting => {
// What is it waiting for? Drop it!
}
CPSComputeStatus::Computing => {
if let Some(waitFor) = waitFor {
if thread.getSettings() == waitFor {
// Wait till thread is finished here.
let (settings, res) = thread.join();
res.and_then(|res| {
if let Err(e) =
Self::storeCPSResultInCache(&file, &settings, &res)
{
self.errors.push(e);
}
None::<()>
});
} else {
// We don't have to wait for this one push
// it back into the stored list of computing
// threads
self.cps_threads.push(thread)
}
} else {
// Put it back in the list of computing threads
self.cps_threads.push(thread)
}
}
CPSComputeStatus::Done(res) => {
if let Some(res) = res
&& let Err(e) = Self::storeCPSResultInCache(&file, &settings, res)
{
self.errors.push(e);
}
}
}
});
}
Err(e) => {
self.errors.push(e);
}
}
}
}