Phono Cartridge Response Measurement Script

[JP]

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.

 

[SeriousSeri]

You are absolutely right. I forgot that I changed the method at some point because of accuracy problems with the end of pilot recognition method! (I was working with spectral energy density) Its been quite some back and forth with those scripts. Pardon me =)

 

[JP]

That’s not encouraging for my hopes and dreams!

 

[SeriousSeri]

That was before i had the right 2-C RIAA correction in place, which made a lot of things a lot more complicated. It should be easier now. We'll see what we can do about it in the future =)

 

[JP]

Got a code snippet for that function?

 

[SeriousSeri]

I actually have around 20 versions of this since i tried a lot of stuff, and i can't remember which one got me closest before i justswitched to sweep start detection.
As i said, this was not tailored for your script, so there were also approaches were i was looking for the start of the pilot tone, but the method for the detection maybe still interesting for you, so i'll drop those snippets aswell, just don't expect too much. Let me know if you want more specific info:

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

 

[JP]

Thank you.

@DrCWO, just converted your scilab script to Python. Not optimized, but working.

 

[DrCWO]

Thank you.

@DrCWO, just converted your scilab script to Python. Not optimized, but working.

Very cool
Is it already included in the measurement script or how did it work?

 

[JP]

Very cool
Is it already included in the measurement script or how did it work?

Not included - just ported it today (I let ChatGPT do most of the work). Going through now to better understand how it works.

 

[DrCWO]

Not included - just ported it today (I let ChatGPT do most of the work). Going through now to better understand how it works.

I wonder if it is easier to build code by myself or try to understand what ChatGPT did

 

[somebodyelse]

Not 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.
Edit: similar to its mathematical workings that look reasonable but have incorrect values if you check them, or URLs that don't actually exist.

 

[JP]

I guess you both missed that I said it works? Of course there were things I had to fix, but it's not complex code so it was a good candidate. I'll happily take the time saver.

 

[JP]

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:

• 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

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 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

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.

 

[DrCWO]

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.

Thank you, this is a big step forward!
Is there already a release I can download an test? As I have to travel maybe I can test it end of next week

 

[JP]

Not yet - still a bit of clean up to do and I need to add the criteria selection for the different test records. Initial testing was just manual tweaking.

 

[JP]

As an update I've promoted version 17.0 to "production" (lol), and I've sent code of the Python version of the file splitter to @DrCWO and @USER that supports TRS-1007 (CA and JVC) and STR-100, and is sample rate agnostic.

 

[JP]

The cutter now works with TRS-1005 50kHz sweeps, though I need to add start-of-sweep detection as the sweep starts 5 seconds after the pilot tone ends. @SeriousSeri's stuff should come in handy for that.

Hopefully we can get this all integrated in to the FR script in January. This will also let me remove the librosa dependency which would be nice as they seem to never support the most current Python release.

 

[USER]

This will also let me remove the librosa dependency which would be nice as they seem to never support the most current Python release.

That's a big deal.

 

[VinyLuke]

Hi All,
I just got the CA-TRS-1007 test record.
How to use it with the script? What settings should I apply in script? I did a first run but I got some nonsence plot with treble falling -21db.
Also should I apply stereo file or delete the silent channel?

EDIT:
ok I reverted to 16.5 from 17.0 and used:
riaamode = 2
riaainv = 1
str100 = 0

Works fine

 

[JP]

What’s the issue with 17? The only change is in the plot layouts, not the data.