- Thread Starter
- #921
Looks like you're looking for 20Hz to show up, rather than looking for the end of the pilot. STR-100 starts at 40Hz and it's way down in the weeds, etc.
def detect_pilot_tone_end(data, sample_rate, pilot_freq=1000, threshold_db=-40):
# Perform STFT
f, t, Zxx = stft(data, fs=sample_rate, nperseg=1024)
# Get the frequency index for the pilot tone
pilot_index = np.argmin(np.abs(f - pilot_freq))
# Convert to dB
Zxx_db = 20 * np.log10(np.abs(Zxx[pilot_index, :]))
# Find the point where the energy at the pilot tone frequency drops below the threshold
end_idx = np.where(Zxx_db < threshold_db)[0]
if end_idx.size > 0:
pilot_tone_end_time = t[end_idx[0]]
else:
pilot_tone_end_time = t[-1] # If not found, use the end of the file
return pilot_tone_end_time
def detect_pilot_tone(channel, samplerate):
# Perform STFT
n_fft = 2048
hop_length = n_fft // 4
D = np.abs(librosa.stft(channel, n_fft=n_fft, hop_length=hop_length))
# Find the frequency index for 1 kHz
freq_bin = int(1000 / (samplerate / float(n_fft)))
pilot_tone_magnitude = D[freq_bin]
threshold = np.max(pilot_tone_magnitude) * 0.5 # 50% of the max magnitude
# Identify time indices where the 1 kHz signal is above threshold
pilot_tone_indices = np.nonzero(pilot_tone_magnitude > threshold)[0]
# Convert time indices to actual time in seconds
pilot_tone_times = librosa.frames_to_time(pilot_tone_indices, sr=samplerate, hop_length=hop_length, n_fft=n_fft)
# Assuming the pilot tone is at the beginning, take the first occurrence
start_time = pilot_tone_times[0] if pilot_tone_times.size > 0 else None
return start_time
def detect_pilot_tone(channel, pilot_tone_freq, samplerate, duration_ms=100):
print("Detecting pilot tone...")
sos = butter(10, [pilot_tone_freq - 10, pilot_tone_freq + 10], btype='bandpass', fs=samplerate, output='sos')
filtered = sosfilt(sos, channel)
envelope = np.abs(hilbert(filtered))
threshold = 0.5 * np.max(envelope)
print("Threshold for pilot tone set to:", threshold)
# Ensure the signal stays above the threshold for the duration of the pilot tone
samples_required = int((duration_ms / 1000.0) * samplerate)
above_threshold = envelope > threshold
continuous_above_threshold = np.where(np.convolve(above_threshold, np.ones(samples_required, dtype=int), 'same') >= samples_required)[0]
if continuous_above_threshold.size > 0:
start_index = continuous_above_threshold[0]
print("Pilot tone detected at index:", start_index)
else:
start_index = None
print("Pilot tone not detected.")
def detect_continuous_tone(channel, freq, samplerate, duration_ms, min_duration=0.5):
sos = butter(10, [freq - 5, freq + 5], btype='bandpass', fs=samplerate, output='sos')
filtered = sosfilt(sos, channel)
envelope = np.abs(hilbert(filtered))
# Compute the threshold as a fraction of the maximum envelope value
threshold = 0.1 * np.max(envelope)
# Convert min_duration to samples
min_duration_samples = int(min_duration * samplerate)
# Find regions where the envelope is above the threshold
above_threshold = envelope > threshold
# Identify regions where the signal is continuously above the threshold
# for at least the minimum duration
groups = np.split(above_threshold, np.where(np.diff(above_threshold) != 0)[0] + 1)
continuous_regions = [group for group in groups if np.all(group) and len(group) > min_duration_samples]
# Find the start of the first continuous region
if continuous_regions:
start_index = np.where(above_threshold == continuous_regions[0][0])[0][0]
else:
start_index = None
def detect_pilot_tone(signal, fs, pilot_tone_freq, pilot_tone_duration, tolerance=0.1):
# Bandpass filter the signal around the pilot tone frequency
filtered = bandpass_filter(signal, pilot_tone_freq - 5, pilot_tone_freq + 5, fs)
# Detect envelope of the filtered signal
envelope = np.abs(hilbert(filtered))
# Thresholding to find where the pilot tone starts and ends
# This can be adjusted based on the specific characteristics of your pilot tone
threshold = np.mean(envelope) * tolerance
above_threshold = envelope > threshold
continuous_tone = np.where(above_threshold)[0]
# Check for a continuous tone of expected duration for the pilot tone
expected_samples = int(pilot_tone_duration * fs)
start_idx = None
end_idx = None
for idx in continuous_tone:
if start_idx is None:
start_idx = idx
elif idx - end_idx > fs * 0.01: # Allowing a gap of 10ms max
# Check if the duration is within the expected range
if end_idx - start_idx > expected_samples * (1 - tolerance) and end_idx - start_idx < expected_samples * (1 + tolerance):
break
else:
start_idx = idx
end_idx = idx
return start_idx, end_idx
def detect_pilot_tone(signal, fs, pilot_tone_freq=1000, min_duration_secs=4):
filtered = bandpass_filter(signal, pilot_tone_freq - 10, pilot_tone_freq + 10, fs, order=2)
envelope = np.abs(hilbert(filtered))
threshold = np.mean(envelope) * 0.1 # Only 10% of the mean to consider fluctuations
continuous_signal = envelope > threshold
# Convert duration to samples
min_duration_samples = min_duration_secs * fs
pilot_tone_start = None
for i in range(len(continuous_signal)):
if continuous_signal[i]:
if pilot_tone_start is None:
pilot_tone_start = i
elif (i - pilot_tone_start) > min_duration_samples:
return pilot_tone_start, i
else:
pilot_tone_start = None
return None, None
Very cool
Very cool
Is it already included in the measurement script or how did it work?
I wonder if it is easier to build code by myself or try to understand what ChatGPT didNot included - just ported it today (I let ChatGPT do most of the work). Going through now to better understand how it works.
And to find out whether ChatGPT made any mistakes... It has a bit of a reputation for spitting out code that looks plausible but doesn't actually do what it's meant to.Not included - just ported it today (I let ChatGPT do most of the work). Going through now to better understand how it works.
Getting the sweeps is not as easy as it seems at first glance. Threshold detections is not sufficient here to get reliable results.The algorithm must also be able to handle speed variations of the turn table and also be robust regarding clicks and pops in the recording. So simple adding N samples to find the end of the sweep (if the Beginning was detected) gives wrong results.
Here is what my script does:
Same can be done for the 20kHz signal to detect the end of the sweep. Here I set the threshold to 5% to get best fit. I tested this approach with 10 needle drops and it works fine so I assume it is quite robust.
- It first I sets up two steep linear phase bandpass filters for 1kHz (16k Samples wide) and 20kHz (1k Samples wide). Linear phase is necessary as filtering should not include any time shift or group delay. Also filters are necessary that don't have much pre- and post-ringing. If someone is interested I can share the impulse responses of the bandbass filters.
View attachment 308293
- Next the scrips reads the recorded WAV from CA-TRS-1007. I assume that the user plays the record from its beginning or drops the needle before the burst followed by the left channel sweep. The recording can be stopped if the third bursts starts playback (see below).
- The left channel (red) gets processed by the following steps:
- Convolve it with the 1kHz Bandpass
- Compute "abs
" so the negative value of the signal gets positive
- Normalize the result to 1 (green). In the graph below I lowered the volume by 12dB to better show both signals.
- Shift the signal half of the filter length (here 8k samples) of the 1kHz Bandpass to compensate for the delay of the linear phase bandpass.
View attachment 308294
Zooming in at the end of the burst it looks like this:
View attachment 308295
We see the green signal has double frequency now. This means at a sample rate of 96kHz the distance between the maxima is 48 samples.
- Next thing is a bit tricky. To get an envelope that can be easily detected I shift the green signal three times by 12 (48/4) samples each and add up the four results. This fills the gaps between the peaks and gives the result below (brown). The end of this signal can be detected easily by thresholding (I use 1% of max) which indicates the beginning of the sweep.
View attachment 308296
It is important that the linear phase bandpass filters have as high attenuation as possible so threshold can be reliably detected. Less attenuation will make thresholding much more difficult. (I tested with other filters befor)
Would be great to see this integrated in the Python Measurement Script to directly read the recorded WAV instead of pre-cut snippets with the sweeps.
An I admit: This can be used for CA-TRS-1007 only and needs RIAA applied before
Edit:
I forgot to mention that I first try to find the beginning of the burst. As the needle drop at the beginning is quite a steep pulse you also get a 1kHz response there. What i did to find the beginning of the burst is to look for 10 consecutive peaks at a distance between 40 and 50 samples so I can be sure I am inside the burst. Starting at this position I look for the end of the burst.
Best DrCWO
Thank you, this is a big step forward!I have it working on variable sample rates and it seems to process flat files just fine, including STR-100. The only hurdle so far is the upper f bandpass to detect the end of the sweep on TRS-1005, unless I go up to 192kHz, otherwise it doesn't detect the end properly as the fundamental goes right up to Fs/2.
As an initial stab I replaced the bandpass filters with 17 order cheby2 IIR as I want the script to be fully self-contained. This mitigated the need to shift the signal and appears to be working fine.
This will also let me remove the librosa dependency which would be nice as they seem to never support the most current Python release.