Python Forum

Hey all,
I have code to generate signal with some peaks and make hilbert - enveloping. But I have some peaks at start in time domain after applying hilbert and then lowpass filter. Any suggestion what I am doing wrong?
Thanks

IMAGE

import numpy as np
import matplotlib.pyplot as plt
import csv
from scipy.signal import hilbert, butter, sosfilt

duration = 1.0  # Duration in seconds
sampling_rate = 30000  # Sampling rate in Hz
num_samples = int(duration * sampling_rate)

time = np.linspace(0, duration, num_samples)
frequency = [25, 134, 234, 95]  # Frequencies of the shaft speed in Hz
amplitude_g_pp = [1.0, 0.8, 0.6, 0.4]  # Peak-to-peak amplitudes of the acceleration in g

acceleration = np.zeros(num_samples)  # Initialize acceleration signal

for freq, amp in zip(frequency, amplitude_g_pp):
    acceleration += amp * np.sin(2 * np.pi * freq * time)

# Add noise at specific frequencies
noise_freqs = [2000,2050,2100,2200,2300, 2350,2400]
noise_amplitude = [0.2, 0.3, 0.5, 0.8, 0.6, 0.4, 0.2]  # Peak-to-peak amplitudes of the acceleration in g

for freq, amp in zip(noise_freqs, noise_amplitude):
    noise_signal = amp * np.sin(2 * np.pi * freq * time)
    acceleration += noise_signal

# Apply Hilbert transform to extract the envelope
envelope = np.abs(hilbert(acceleration))

# Apply bandpass filter to the envelope signal
lowcut = 500
highcut = 14900
nyquist = 0.5 * sampling_rate
low = lowcut / nyquist
high = highcut / nyquist
sos = butter(4, [low, high], btype='band', output='sos')
filtered_envelope = sosfilt(sos, envelope)

# Save envelope values to a CSV file
csv_file = "envelope_data.csv"
header = ["Envelope (g)"]

with open(csv_file, 'w', newline='') as file:
    writer = csv.writer(file)
    writer.writerows(np.around(filtered_envelope, decimals=8).reshape(-1, 1))

print("Envelope values saved to", csv_file)

# Plotting time-domain signal
plt.subplot(3, 2, 1)
plt.plot(time, acceleration)
plt.xlabel('Time (s)')
plt.ylabel('Acceleration (g)')
plt.title('Acceleration Signal')
plt.grid(True)

# Plotting FFT spectrum of the acceleration signal
fft_acc = np.fft.fft(acceleration)
freq_acc = np.fft.fftfreq(num_samples, 1 / sampling_rate)
amplitude_acc = np.abs(fft_acc) / num_samples * 2  # Normalize the FFT

# Select the desired frequency range for FFT plot of acceleration signal
fft_acc_freq_range = (freq_acc >= 10) & (freq_acc <= 15000)
freq_acc = freq_acc[fft_acc_freq_range]
amplitude_acc = amplitude_acc[fft_acc_freq_range]

plt.subplot(3, 2, 2)
plt.plot(freq_acc, amplitude_acc)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude (g)')
plt.title('FFT of Acceleration Signal (10-15000 Hz)')
plt.grid(True)

# Plotting the envelope of the acceleration signal
plt.subplot(3, 2, 3)
plt.plot(time, filtered_envelope)
plt.xlabel('Time (s)')
plt.ylabel('Envelope (g)')
plt.title('Envelope of Acceleration Signal')
plt.grid(True)

# Plotting the contour of the upper wave of the filtered envelope signal
upper_wave_contour = np.maximum(filtered_envelope, 0)

# Apply FFT to the windowed contour signal
fft_contour = np.fft.fft(upper_wave_contour)

# Apply lowpass filter in frequency domain
cutoff_freq = 1000
freq = np.fft.fftfreq(len(fft_contour), 1 / sampling_rate)
fft_contour[np.abs(freq) > cutoff_freq] = 0

# Obtain the filtered contour by inverse FFT
filtered_contour = np.fft.ifft(fft_contour).real

# Plotting the filtered contour of the upper wave of the filtered envelope signal
plt.subplot(3, 2, 4)
plt.plot(time, filtered_contour)
plt.xlabel('Time (s)')
plt.ylabel('Contour')
plt.title('Envelope Signal - lowpass')
plt.grid(True)

# Plotting FFT spectrum of the filtered envelope signal
fft_env = np.fft.fft(filtered_contour)
freq_env = np.fft.fftfreq(num_samples, 1 / sampling_rate)
amplitude_env = np.abs(fft_env) / num_samples * 2  # Normalize the FFT

# Select the desired frequency range for FFT plot of the filtered envelope signal
fft_env_freq_range = (freq_env >= 10) & (freq_env <= 1000)
freq_env = freq_env[fft_env_freq_range]
amplitude_env = amplitude_env[fft_env_freq_range]

plt.subplot(3, 2, 5)
plt.plot(freq_env, amplitude_env)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude (g)')
plt.title('FFT of Filtered Envelope Signal (10-1000 Hz)')
plt.grid(True)

# Calculate velocity from acceleration using cumulative integration (cumtrapz)
velocity = np.cumsum(acceleration) / sampling_rate * 1000  * 9.8   # Convert to mm/s

# Plotting FFT spectrum of the velocity signal
fft_vel = np.fft.fft(velocity)
freq_vel = np.fft.fftfreq(num_samples, 1 / sampling_rate)
amplitude_vel = np.abs(fft_vel) / num_samples * 2  # Normalize the FFT

# Select the desired frequency range for FFT plot of the velocity signal
fft_vel_freq_range = (freq_vel >= 10) & (freq_vel <= 1000)
freq_vel = freq_vel[fft_vel_freq_range]
amplitude_vel = amplitude_vel[fft_vel_freq_range]

plt.subplot(3, 2, 6)
plt.plot(freq_vel, amplitude_vel)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude (mm/s)')
plt.title('FFT of Velocity Signal (10-1000 Hz)')
plt.grid(True)

# Adjusting subplots layout
plt.tight_layout()

# Display the plot
plt.show()