(band, timeSignal, samplingRate, numSamples,
numChannels, timeLength, suppressWarnings=True)
| 1247 | |
| 1248 | # @njit |
| 1249 | def _Lundeby_correction(band, timeSignal, samplingRate, numSamples, |
| 1250 | numChannels, timeLength, suppressWarnings=True): |
| 1251 | returnTuple = (np.float32(0), np.float32(0), np.int32(0), np.float32(0)) |
| 1252 | timeSignal, sampleShift = _circular_time_shift(timeSignal) |
| 1253 | if sampleShift is None: |
| 1254 | return returnTuple |
| 1255 | |
| 1256 | numSamples -= sampleShift # discount shifted samples |
| 1257 | numParts = 5 # number of parts per 10 dB decay. N = any([3, 10]) |
| 1258 | dBtoNoise = 7 # stop point 10 dB above first estimated background noise |
| 1259 | useDynRange = 15 # dynamic range |
| 1260 | |
| 1261 | # Window length - 10 to 50 ms, longer periods for lower frequencies and vice versa |
| 1262 | repeat = True |
| 1263 | i = 0 |
| 1264 | winTimeLength = 0.01 |
| 1265 | while repeat: # loop to find proper winTimeLength |
| 1266 | winTimeLength = winTimeLength + 0.01*i |
| 1267 | # 1) local time average: |
| 1268 | blockSamples = int(winTimeLength * samplingRate) |
| 1269 | timeWinData, timeVecWin = _level_profile(timeSignal, samplingRate, |
| 1270 | numSamples, numChannels, |
| 1271 | blockSamples) |
| 1272 | |
| 1273 | # 2) estimate noise from h^2_averaged(block): |
| 1274 | bgNoiseLevel = 10 * \ |
| 1275 | np.log10( |
| 1276 | np.mean(timeWinData[-int(timeWinData.size/10):])) |
| 1277 | |
| 1278 | # 3) Calculate preliminar slope |
| 1279 | startIdx = np.argmax(np.abs(timeWinData/np.max(np.abs(timeWinData)))) |
| 1280 | stopIdx = startIdx + np.where(10*np.log10(timeWinData[startIdx+1:]) |
| 1281 | >= bgNoiseLevel + dBtoNoise)[0][-1] |
| 1282 | dynRange = 10*np.log10(timeWinData[stopIdx]) \ |
| 1283 | - 10*np.log10(timeWinData[startIdx]) |
| 1284 | if (stopIdx == startIdx) or (dynRange > -5)[0]: |
| 1285 | if not suppressWarnings: |
| 1286 | print(band, "[Hz] band: SNR too low for the preliminar slope", |
| 1287 | "calculation.") |
| 1288 | # return returnTuple |
| 1289 | |
| 1290 | # X*c = EDC (energy decaying curve) |
| 1291 | X = np.ones((stopIdx-startIdx, 2), dtype=np.float32) |
| 1292 | X[:, 1] = timeVecWin[startIdx:stopIdx, 0] |
| 1293 | c = np.linalg.lstsq(X, 10*np.log10(timeWinData[startIdx:stopIdx]), |
| 1294 | rcond=-1)[0] |
| 1295 | |
| 1296 | if (c[1] == 0)[0] or np.isnan(c).any(): |
| 1297 | if not suppressWarnings: |
| 1298 | print(band, "[Hz] band: regression failed. T would be inf.") |
| 1299 | # return returnTuple |
| 1300 | |
| 1301 | # 4) preliminary intersection |
| 1302 | crossingPoint = (bgNoiseLevel - c[0]) / c[1] # [s] |
| 1303 | if (crossingPoint > 2*(timeLength + sampleShift/samplingRate))[0]: |
| 1304 | if not suppressWarnings: |
| 1305 | print(band, "[Hz] band: preliminary intersection point between", |
| 1306 | "bgNoiseLevel and the decay slope greater than signal length.") |
no test coverage detected