(原创)超声波/雷达信号处理仿真与FPGA加速
一、标准信号产生与接收
1、标准信号的产生与接收
1、要求首先使用10Mhz Sine波进行发送一次持续10个Sine周期的脉冲信号模拟超声波信号发送。在不发送的时候，就保持0。
2、假设物体距离探头的距离为10米，当前声速为1000m/s，请计算超声波碰到物体后弹回来的总时间作为延时。然后在发射的信号前增加这个延时作为模拟接收到的原始信号。
3、将发射和接受到的信号以图像的形式进行显示出来。
4、计算接收到的信号的FFT，频率范围为0到5倍信号原始频率(10Mhz)。
5、对接收到的信号与发射的信号脉冲进行互相关计算，然后找到最大值，计算发射和接收之间的延时。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate
# Constants
f = 10e6 # Frequency of the sine wave (10 MHz)
distance = 10 # Distance to the object in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Time settings
pulse_duration = 10 / f # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate transmitted signal (10 MHz sine wave)
transmitted_signal = np.zeros_like(t)
transmitted_signal[:pulse_samples] = np.sin(2 * np.pi * f * t[:pulse_samples])
# Calculate round trip time
round_trip_time = 2 * distance / speed_of_sound # Total time for the sound to travel to the object and back
# Simulate received signal with delay
received_signal = np.zeros_like(t)
received_delay_index = int(round_trip_time / T)
received_signal[received_delay_index:received_delay_index + pulse_samples] = transmitted_signal[:pulse_samples]
# Plot transmitted and received signals
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.title('Transmitted and Received Signals')
plt.plot(t, transmitted_signal, label='Transmitted Signal', color='blue')
plt.plot(t, received_signal, label='Received Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
# FFT of received signal
N = len(received_signal)
frequencies = fftfreq(N, T)[:N // 2]
received_signal_fft = fft(received_signal)
# Plot FFT
plt.subplot(2, 1, 2)
plt.title('FFT of Received Signal')
plt.plot(frequencies, 2.0/N * np.abs(received_signal_fft[:N // 2]), color='green')
plt.xlim(0, 5 * f) # Limit to 5 times the original frequency
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()
# Cross-correlation
correlation = correlate(received_signal, transmitted_signal, mode='full')
lag = np.arange(-len(transmitted_signal) + 1, len(received_signal))
# Find the delay
max_corr_index = np.argmax(correlation)
max_delay = lag[max_corr_index] * T
print(f'Maximum correlation delay: {max_delay:.6f} seconds')运行结果：
二、添加干扰和直流偏置
1、发射信号的幅值为10V。
2、模拟的接收信号，添加随机噪声，以尽可能模拟真实情况。然后模拟的接收信号全部添加100V的直流偏置。
3、画出互相关的图像。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate
# Constants
f = 10e6 # Frequency of the sine wave (10 MHz)
distance = 10 # Distance to the object in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Time settings
pulse_duration = 10 / f # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate transmitted signal (10 MHz sine wave)
transmitted_signal = np.zeros_like(t)
transmitted_signal[:pulse_samples] = 10 * np.sin(2 * np.pi * f * t[:pulse_samples]) # Amplitude 10V
# Calculate round trip time
round_trip_time = 2 * distance / speed_of_sound # Total time for the sound to travel to the object and back
# Simulate received signal with delay
received_signal = np.zeros_like(t)
received_delay_index = int(round_trip_time / T)
received_signal[received_delay_index:received_delay_index + pulse_samples] = transmitted_signal[:pulse_samples]
# Add random noise to the received signal
noise = np.random.normal(0, 1, received_signal.shape) # Gaussian noise
received_signal += noise
# Add DC offset of 100V
received_signal += 100
# Plot transmitted and received signals
plt.figure(figsize=(12, 6))
plt.subplot(3, 1, 1)
plt.title('Transmitted and Received Signals')
plt.plot(t, transmitted_signal, label='Transmitted Signal', color='blue')
plt.plot(t, received_signal, label='Received Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
# FFT of received signal
N = len(received_signal)
frequencies = fftfreq(N, T)[:N // 2]
received_signal_fft = fft(received_signal)
# Plot FFT
plt.subplot(3, 1, 2)
plt.title('FFT of Received Signal')
plt.plot(frequencies, 2.0/N * np.abs(received_signal_fft[:N // 2]), color='green')
plt.xlim(0, 5 * f) # Limit to 5 times the original frequency
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
# Cross-correlation
correlation = correlate(received_signal, transmitted_signal, mode='full')
lag = np.arange(-len(transmitted_signal) + 1, len(received_signal))
# Plot cross-correlation
plt.subplot(3, 1, 3)
plt.title('Cross-Correlation')
plt.plot(lag * T, correlation, color='purple')
plt.xlabel('Lag (s)')
plt.ylabel('Correlation')
plt.tight_layout()
plt.show()
# Find the delay
max_corr_index = np.argmax(correlation)
max_delay = lag[max_corr_index] * T
print(f'Maximum correlation delay: {max_delay:.6f} seconds')运行结果：
三、滤波
1、接收到信号之后，使用切比雪夫I滤波器实现一个高通滤波器，滤除直流分量。
2、高通滤波器之后，进行低通滤波器，已知信号频率为10Mhz，设计合适的滤波器。
3、对滤波前和滤波后的接收信号进行显示，同时绘制FFT图像。
4、对滤波后的信号进行互相关计算和图像显示。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate, cheby1, lfilter
# Constants
f = 10e6 # Frequency of the sine wave (10 MHz)
distance = 10 # Distance to the object in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Time settings
pulse_duration = 10 / f # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate transmitted signal (10 MHz sine wave)
transmitted_signal = np.zeros_like(t)
transmitted_signal[:pulse_samples] = 10 * np.sin(2 * np.pi * f * t[:pulse_samples]) # Amplitude 10V
# Calculate round trip time
round_trip_time = 2 * distance / speed_of_sound # Total time for the sound to travel to the object and back
# Simulate received signal with delay
received_signal = np.zeros_like(t)
received_delay_index = int(round_trip_time / T)
received_signal[received_delay_index:received_delay_index + pulse_samples] = transmitted_signal[:pulse_samples]
# Add random noise to the received signal
noise = np.random.normal(0, 1, received_signal.shape) # Gaussian noise
received_signal += noise
# Add DC offset of 100V
received_signal += 100
# Design Chebyshev I high-pass filter to remove DC component
nyquist = 0.5 * sampling_rate
cutoff_high = 1e5 # High-pass cutoff frequency (100 kHz)
order = 1
b, a = cheby1(order, 0.5, cutoff_high / nyquist, btype='highpass')
filtered_signal_high = lfilter(b, a, received_signal)
# Design Chebyshev I low-pass filter to remove high-frequency noise
cutoff_low = 15e6 # Low-pass cutoff frequency (15 MHz)
b_low, a_low = cheby1(order, 0.5, cutoff_low / nyquist, btype='low')
filtered_signal_low = lfilter(b_low, a_low, filtered_signal_high)
# Plot received and filtered signals
plt.figure(figsize=(12, 8))
plt.subplot(4, 1, 1)
plt.title('Received Signal with Noise and DC Offset')
plt.plot(t, received_signal, label='Received Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 2)
plt.title('Filtered Signal (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low, label='Filtered Signal', color='green')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
# FFT of filtered signal
N = len(filtered_signal_low)
frequencies = fftfreq(N, T)[:N // 2]
filtered_signal_fft = fft(filtered_signal_low)
# Plot FFT of filtered signal
plt.subplot(4, 1, 3)
plt.title('FFT of Filtered Signal')
plt.plot(frequencies, 2.0/N * np.abs(filtered_signal_fft[:N // 2]), color='blue')
plt.xlim(0, 5 * f) # Limit to 5 times the original frequency
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
# Cross-correlation
correlation = correlate(filtered_signal_low, transmitted_signal, mode='full')
lag = np.arange(-len(transmitted_signal) + 1, len(filtered_signal_low))
# Plot cross-correlation
plt.subplot(4, 1, 4)
plt.title('Cross-Correlation of Filtered Signal')
plt.plot(lag * T, correlation, color='purple')
plt.xlabel('Lag (s)')
plt.ylabel('Correlation')
plt.tight_layout()
plt.show()
# Find the delay
max_corr_index = np.argmax(correlation)
max_delay = lag[max_corr_index] * T
print(f'Maximum correlation delay: {max_delay:.6f} seconds')
截止到目前，基本都功能已经实现了，但是我们仍然可以进行改进。
目前已知问题：
1、滤波后的波形在初始位置有个瞬态响应。
2、按理说滤波器会改变高频时的相位，滤波器的设计需要小心。
3、没有考虑到多普勒效应。
4、多目标的距离检测没有实现。
四、多点探测
在多点检测时，可能会发生混叠现象，解决方法是：1、更换脉冲形状，改为不定频脉冲，例如LFM脉冲。2、降低脉宽，但可能会带来信噪比的问题。
1、首先修改脉冲为LFM。
2、增加一个物体的模拟，即原先距离10米的物体不变，增加一个距离50米的物体。其他条件保持不变。
3、能够通过互相关识别到多个物体(超过最大峰值50%的峰均认为是有效的物体距离)。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate, cheby1, lfilter
# Constants
f0 = 10e6 # Start frequency of the LFM pulse (10 MHz)
distance1 = 10 # Distance to the first object in meters
distance2 = 50 # Distance to the second object in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Time settings
pulse_duration = 10 / f0 # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate LFM pulse
k = (1e10) # Chirp rate
lfm_pulse = np.zeros_like(t)
lfm_pulse[:pulse_samples] = np.sin(2 * np.pi * (f0 * t[:pulse_samples] + 0.5 * k * t[:pulse_samples]**2)) # LFM pulse
# Simulate received signals with delays
received_signal = np.zeros_like(t)
received_delay_index1 = int((2 * distance1) / speed_of_sound / T)
received_delay_index2 = int((2 * distance2) / speed_of_sound / T)
received_signal[received_delay_index1:received_delay_index1 + pulse_samples] = lfm_pulse[:pulse_samples]
received_signal[received_delay_index2:received_delay_index2 + pulse_samples] += lfm_pulse[:pulse_samples]
# Add random noise to the received signal (lower standard deviation for realistic simulation)
noise_std_dev = 0.1 # Standard deviation of noise
noise = np.random.normal(0, noise_std_dev, received_signal.shape) # Gaussian noise
received_signal += noise
# Add DC offset of 100V
received_signal += 100
# Design Chebyshev I high-pass filter to remove DC component
nyquist = 0.5 * sampling_rate
cutoff_high = 1e5 # High-pass cutoff frequency (100 kHz)
order = 1
b, a = cheby1(order, 0.5, cutoff_high / nyquist, btype='highpass')
filtered_signal_high = lfilter(b, a, received_signal)
# Design Chebyshev I low-pass filter to remove high-frequency noise
cutoff_low = 15e6 # Low-pass cutoff frequency (15 MHz)
b_low, a_low = cheby1(order, 0.5, cutoff_low / nyquist, btype='low')
filtered_signal_low = lfilter(b_low, a_low, filtered_signal_high)
# Plot received and filtered signals
plt.figure(figsize=(12, 8))
plt.subplot(4, 1, 1)
plt.title('Received Signal with Noise and DC Offset')
plt.plot(t, received_signal, label='Received Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 2)
plt.title('Filtered Signal (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low, label='Filtered Signal', color='green')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
# FFT of filtered signal
N = len(filtered_signal_low)
frequencies = fftfreq(N, T)[:N // 2]
filtered_signal_fft = fft(filtered_signal_low)
# Plot FFT of filtered signal
plt.subplot(4, 1, 3)
plt.title('FFT of Filtered Signal')
plt.plot(frequencies, 2.0/N * np.abs(filtered_signal_fft[:N // 2]), color='blue')
plt.xlim(0, 5 * f0) # Limit to 5 times the original frequency
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
# Cross-correlation
correlation = correlate(filtered_signal_low, lfm_pulse, mode='full')
lag = np.arange(-len(lfm_pulse) + 1, len(filtered_signal_low))
# Plot cross-correlation
plt.subplot(4, 1, 4)
plt.title('Cross-Correlation of Filtered Signal')
plt.plot(lag * T, correlation, color='purple')
plt.xlabel('Lag (s)')
plt.ylabel('Correlation')
plt.tight_layout()
plt.show()
# Identify peaks in cross-correlation
threshold = 0.5 * np.max(correlation) # 50% of the maximum peak
peaks = np.where(correlation > threshold)[0]
# Calculate distances for detected objects
detected_distances = []
for peak in peaks:
distance_detected = (peak - len(lfm_pulse) + 1) * T * speed_of_sound / 2
detected_distances.append(distance_detected)
print(f'Detected distances: {detected_distances}')运行结果：
但是算出来了很多点：
PS C:\Users\ChanRa1n> & C:/Users/ChanRa1n/AppData/Local/Programs/Python/Python37/python.exe c:/Users/ChanRa1n/Desktop/ultrasonic.py Detected distances: [-9.999999999999999e-05, -4.9999999999999996e-05, 0.0, 5e-06, 1e-05, 1.5000000000000002e-05, 2e-05, 2.4999999999999998e-05, 3.0000000000000004e-05, 3.5000000000000004e-05, 4e-05, 4.4999999999999996e-05, 4.9999999999999996e-05, 5.5e-05, 6.000000000000001e-05, 6.5e-05, 7.000000000000001e-05, 7.5e-05, 8e-05, 8.5e-05, 8.999999999999999e-05, 9.5e-05, 9.999999999999999e-05, 0.000105, 0.00011, 0.000115, 0.00012000000000000002, 0.000125, 0.00013, 0.000135, 0.00014000000000000001, 0.00014500000000000003, 0.00015, 0.000155, 0.00016, 0.000165, 0.00017, 0.000175, 0.00017999999999999998, 0.000185, 0.00019, 0.00019500000000000002, 0.00019999999999999998, 0.000205, 0.00021, 0.000215, 0.00022, 0.00022500000000000002, 0.00023, 0.000235, 0.00024000000000000003, 0.000245, 0.00025, 0.000255, 0.00026, 0.000265, 0.00027, 0.000275, 0.00028000000000000003, 0.00028500000000000004, 0.00029000000000000006, 0.00029499999999999996, 0.0003, 0.000305, 0.00031, 0.000315, 0.00032, 0.000325, 0.00033, 0.000335, 0.00034, 0.00034500000000000004, 0.00035, 0.00037, 0.000375, 0.00038, 0.00038500000000000003, 0.00039000000000000005, 0.000395, 0.00039999999999999996, 0.00042500000000000003, 0.00043, 0.000435, 0.00044, 0.00044500000000000003, 0.000485, 0.00049, 0.000495, 9.99975, 9.999755, 9.9998, 9.999805, 9.999849999999999, 9.999855, 9.999895, 9.9999, 9.999905, 9.999945, 9.99995, 9.999955, 9.999995, 10.0, 10.000005000000002, 10.00001, 10.000045, 10.00005, 10.000055000000001, 10.000095, 10.0001, 10.000105, 10.000150000000001, 10.000155, 10.000200000000001, 10.000205, 49.9998, 49.999805, 49.999849999999995, 49.999855000000004, 49.999895, 49.9999, 49.999905000000005, 49.999945000000004, 49.99995, 49.999955, 49.999995, 50.0, 50.000005, 50.00001, 50.000045, 50.00005, 50.000055, 50.000099999999996, 50.000105000000005, 50.00015, 50.00015500000001, 50.0002, 50.000205]
但是很明显存在问题，所以我修改了代码，首先找到所有极大值，然后在极大值中进行设定阈值，超过阈值的部分方可被认为是有效的。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate, cheby1, lfilter, find_peaks
# Constants
f0 = 10e6 # Start frequency of the LFM pulse (10 MHz)
distance1 = 10 # Distance to the first object in meters
distance2 = 50 # Distance to the second object in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Time settings
pulse_duration = 10 / f0 # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate LFM pulse
k = (1e10) # Chirp rate
lfm_pulse = np.zeros_like(t)
lfm_pulse[:pulse_samples] = np.sin(2 * np.pi * (f0 * t[:pulse_samples] + 0.5 * k * t[:pulse_samples]**2)) # LFM pulse
# Simulate received signals with delays
received_signal = np.zeros_like(t)
received_delay_index1 = int((2 * distance1) / speed_of_sound / T)
received_delay_index2 = int((2 * distance2) / speed_of_sound / T)
received_signal[received_delay_index1:received_delay_index1 + pulse_samples] = lfm_pulse[:pulse_samples]
received_signal[received_delay_index2:received_delay_index2 + pulse_samples] += lfm_pulse[:pulse_samples]
# Add random noise to the received signal (lower standard deviation for realistic simulation)
noise_std_dev = 0.1 # Standard deviation of noise
noise = np.random.normal(0, noise_std_dev, received_signal.shape) # Gaussian noise
received_signal += noise
# Add DC offset of 100V
received_signal += 100
# Design Chebyshev I high-pass filter to remove DC component
nyquist = 0.5 * sampling_rate
cutoff_high = 1e5 # High-pass cutoff frequency (100 kHz)
order = 1
b, a = cheby1(order, 0.5, cutoff_high / nyquist, btype='highpass')
filtered_signal_high = lfilter(b, a, received_signal)
# Design Chebyshev I low-pass filter to remove high-frequency noise
cutoff_low = 15e6 # Low-pass cutoff frequency (15 MHz)
b_low, a_low = cheby1(order, 0.5, cutoff_low / nyquist, btype='low')
filtered_signal_low = lfilter(b_low, a_low, filtered_signal_high)
# Plot received and filtered signals
plt.figure(figsize=(12, 8))
plt.subplot(4, 1, 1)
plt.title('Received Signal with Noise and DC Offset')
plt.plot(t, received_signal, label='Received Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 2)
plt.title('Filtered Signal (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low, label='Filtered Signal', color='green')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
# FFT of filtered signal
N = len(filtered_signal_low)
frequencies = fftfreq(N, T)[:N // 2]
filtered_signal_fft = fft(filtered_signal_low)
# Plot FFT of filtered signal
plt.subplot(4, 1, 3)
plt.title('FFT of Filtered Signal')
plt.plot(frequencies, 2.0/N * np.abs(filtered_signal_fft[:N // 2]), color='blue')
plt.xlim(0, 5 * f0) # Limit to 5 times the original frequency
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
# Cross-correlation
correlation = correlate(filtered_signal_low, lfm_pulse, mode='full')
lag = np.arange(-len(lfm_pulse) + 1, len(filtered_signal_low))
# Plot cross-correlation
plt.subplot(4, 1, 4)
plt.title('Cross-Correlation of Filtered Signal')
plt.plot(lag * T, correlation, color='purple')
plt.xlabel('Lag (s)')
plt.ylabel('Correlation')
plt.tight_layout()
plt.show()
# Identify peaks in cross-correlation
peaks, _ = find_peaks(correlation) # Find peaks in the correlation
threshold = 0.6 * np.max(correlation) # 60% of the maximum peak
valid_peaks = peaks[correlation[peaks] > threshold] # Filter peaks based on threshold
# Calculate distances for detected objects
detected_distances = []
for peak in valid_peaks:
distance_detected = (peak - len(lfm_pulse) + 1) * T * speed_of_sound / 2
detected_distances.append(distance_detected)
print(f'Detected distances: {detected_distances}')这个时候结果好了一点：
Detected distances: [-4.9999999999999996e-05, 0.0, 1e-05, 6.5e-05, 9.999849999999999, 9.9999, 9.99995, 10.0, 10.00005, 10.0001, 10.000150000000001, 49.999849999999995, 49.9999, 49.99995, 50.0, 50.00005, 50.000099999999996, 50.00015, 50.0002]
进一步优化：
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate, cheby1, lfilter, find_peaks
# Constants
f0 = 10e6 # Start frequency of the LFM pulse (10 MHz)
distance1 = 10 # Distance to the first object in meters
distance2 = 50 # Distance to the second object in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Time settings
pulse_duration = 10 / f0 # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate LFM pulse
k = (1e10) # Chirp rate
lfm_pulse = np.zeros_like(t)
lfm_pulse[:pulse_samples] = np.sin(2 * np.pi * (f0 * t[:pulse_samples] + 0.5 * k * t[:pulse_samples]**2)) # LFM pulse
# Simulate received signals with delays
received_signal = np.zeros_like(t)
received_delay_index1 = int((2 * distance1) / speed_of_sound / T)
received_delay_index2 = int((2 * distance2) / speed_of_sound / T)
received_signal[received_delay_index1:received_delay_index1 + pulse_samples] = lfm_pulse[:pulse_samples]
received_signal[received_delay_index2:received_delay_index2 + pulse_samples] += lfm_pulse[:pulse_samples]
# Add random noise to the received signal (lower standard deviation for realistic simulation)
noise_std_dev = 0.1 # Standard deviation of noise
noise = np.random.normal(0, noise_std_dev, received_signal.shape) # Gaussian noise
received_signal += noise
# Add DC offset of 100V
received_signal += 100
# Design Chebyshev I high-pass filter to remove DC component
nyquist = 0.5 * sampling_rate
cutoff_high = 1e5 # High-pass cutoff frequency (100 kHz)
order = 1
b, a = cheby1(order, 0.5, cutoff_high / nyquist, btype='highpass')
filtered_signal_high = lfilter(b, a, received_signal)
# Design Chebyshev I low-pass filter to remove high-frequency noise
cutoff_low = 15e6 # Low-pass cutoff frequency (15 MHz)
b_low, a_low = cheby1(order, 0.5, cutoff_low / nyquist, btype='low')
filtered_signal_low = lfilter(b_low, a_low, filtered_signal_high)
# Plot received and filtered signals
plt.figure(figsize=(12, 8))
plt.subplot(4, 1, 1)
plt.title('Received Signal with Noise and DC Offset')
plt.plot(t, received_signal, label='Received Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 2)
plt.title('Filtered Signal (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low, label='Filtered Signal', color='green')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
# FFT of filtered signal
N = len(filtered_signal_low)
frequencies = fftfreq(N, T)[:N // 2]
filtered_signal_fft = fft(filtered_signal_low)
# Plot FFT of filtered signal
plt.subplot(4, 1, 3)
plt.title('FFT of Filtered Signal')
plt.plot(frequencies, 2.0/N * np.abs(filtered_signal_fft[:N // 2]), color='blue')
plt.xlim(0, 5 * f0) # Limit to 5 times the original frequency
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
# Cross-correlation
correlation = correlate(filtered_signal_low, lfm_pulse, mode='full')
lag = np.arange(-len(lfm_pulse) + 1, len(filtered_signal_low))
# Plot cross-correlation
plt.subplot(4, 1, 4)
plt.title('Cross-Correlation of Filtered Signal')
plt.plot(lag * T, correlation, color='purple')
plt.xlabel('Lag (s)')
plt.ylabel('Correlation')
plt.tight_layout()
plt.show()
# Identify peaks in cross-correlation
peaks, _ = find_peaks(correlation) # Find peaks in the correlation
threshold = 0.8 * np.max(correlation) # 80% of the maximum peak
valid_peaks = peaks[correlation[peaks] > threshold] # Filter peaks based on threshold
# Calculate distances for detected objects
detected_distances = []
tolerance = 0.01 # Tolerance for detecting unique distances
last_distance = None
for peak in valid_peaks:
distance_detected = abs((peak - len(lfm_pulse) + 1) * T * speed_of_sound / 2)
if last_distance is None or abs(distance_detected - last_distance) > tolerance:
detected_distances.append(distance_detected)
last_distance = distance_detected
print(f'Detected distances: {detected_distances}')此时的结果：
PS C:\Users\ChanRa1n> & C:/Users/ChanRa1n/AppData/Local/Programs/Python/Python37/python.exe c:/Users/ChanRa1n/Desktop/ultrasonic.py Detected distances: [9.9999, 49.9999]
和我们设计的10米和50米非常接近了。为了减少FFT计算的频谱泄露问题，添加汉明窗函数。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate, cheby1, lfilter, find_peaks, hamming
# Constants
f0 = 10e6 # Start frequency of the LFM pulse (10 MHz)
distance1 = 10 # Distance to the first object in meters
distance2 = 50 # Distance to the second object in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Time settings
pulse_duration = 10 / f0 # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate LFM pulse
k = (1e10) # Chirp rate
lfm_pulse = np.zeros_like(t)
lfm_pulse[:pulse_samples] = np.sin(2 * np.pi * (f0 * t[:pulse_samples] + 0.5 * k * t[:pulse_samples]**2)) # LFM pulse
# Simulate received signals with delays
received_signal = np.zeros_like(t)
received_delay_index1 = int((2 * distance1) / speed_of_sound / T)
received_delay_index2 = int((2 * distance2) / speed_of_sound / T)
received_signal[received_delay_index1:received_delay_index1 + pulse_samples] = lfm_pulse[:pulse_samples]
received_signal[received_delay_index2:received_delay_index2 + pulse_samples] += lfm_pulse[:pulse_samples]
# Add random noise to the received signal (lower standard deviation for realistic simulation)
noise_std_dev = 0.1 # Standard deviation of noise
noise = np.random.normal(0, noise_std_dev, received_signal.shape) # Gaussian noise
received_signal += noise
# Add DC offset of 100V
received_signal += 100
# Design Chebyshev I high-pass filter to remove DC component
nyquist = 0.5 * sampling_rate
cutoff_high = 1e5 # High-pass cutoff frequency (100 kHz)
order = 1
b, a = cheby1(order, 0.5, cutoff_high / nyquist, btype='highpass')
filtered_signal_high = lfilter(b, a, received_signal)
# Design Chebyshev I low-pass filter to remove high-frequency noise
cutoff_low = 15e6 # Low-pass cutoff frequency (15 MHz)
b_low, a_low = cheby1(order, 0.5, cutoff_low / nyquist, btype='low')
filtered_signal_low = lfilter(b_low, a_low, filtered_signal_high)
# Plot received and filtered signals
plt.figure(figsize=(12, 8))
plt.subplot(4, 1, 1)
plt.title('Received Signal with Noise and DC Offset')
plt.plot(t, received_signal, label='Received Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 2)
plt.title('Filtered Signal (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low, label='Filtered Signal', color='green')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
# Apply Hamming window to the filtered signal before FFT
window = hamming(len(filtered_signal_low))
windowed_signal = filtered_signal_low * window
# FFT of windowed signal
N = len(windowed_signal)
frequencies = fftfreq(N, T)[:N // 2]
filtered_signal_fft = fft(windowed_signal)
# Plot FFT of filtered signal
plt.subplot(4, 1, 3)
plt.title('FFT of Filtered Signal with Hamming Window')
plt.plot(frequencies, 2.0/N * np.abs(filtered_signal_fft[:N // 2]), color='blue')
plt.xlim(0, 5 * f0) # Limit to 5 times the original frequency
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
# Cross-correlation
correlation = correlate(filtered_signal_low, lfm_pulse, mode='full')
lag = np.arange(-len(lfm_pulse) + 1, len(filtered_signal_low))
# Plot cross-correlation
plt.subplot(4, 1, 4)
plt.title('Cross-Correlation of Filtered Signal')
plt.plot(lag * T, correlation, color='purple')
plt.xlabel('Lag (s)')
plt.ylabel('Correlation')
plt.tight_layout()
plt.show()
# Identify peaks in cross-correlation
peaks, _ = find_peaks(correlation) # Find peaks in the correlation
threshold = 0.8 * np.max(correlation) # 80% of the maximum peak
valid_peaks = peaks[correlation[peaks] > threshold] # Filter peaks based on threshold
# Calculate distances for detected objects
detected_distances = []
tolerance = 0.01 # Tolerance for detecting unique distances
last_distance = None
for peak in valid_peaks:
distance_detected = abs((peak - len(lfm_pulse) + 1) * T * speed_of_sound / 2)
if last_distance is None or abs(distance_detected - last_distance) > tolerance:
detected_distances.append(distance_detected)
last_distance = distance_detected
print(f'Detected distances: {detected_distances}')运行结果：
PS C:\Users\ChanRa1n> & C:/Users/ChanRa1n/AppData/Local/Programs/Python/Python37/python.exe c:/Users/ChanRa1n/Desktop/ultrasonic.py Detected distances: [9.9999, 49.99995]
之所以会有误差，主要是来自两方面：1、在计算过程中，浮点数的精度可能导致结果略微偏离理论值。2、滤波器会稍微改变信号的相位。
由于误差在可接受范围内，故此处不进行修正(后面要在FPGA上实现，高精度会大幅度增加面积)。
五、多阵元单目标探测
假设存在第二个阵元，与第一个阵元的直线距离为1米，在本次实验中始终仅进行接收。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate, cheby1, lfilter
# Constants
f = 10e6 # Frequency of the sine wave (10 MHz)
distance_vertical = 10 # Vertical distance to the object in meters
distance_horizontal = 1 # Horizontal distance between array 1 and array 2 in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Calculate the actual distances to the object for each array
distance_1 = distance_vertical # Distance from array 1 to the object
distance_2 = np.sqrt(distance_vertical**2 + distance_horizontal**2) # Distance from array 2 to the object
# Time settings
pulse_duration = 10 / f # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate transmitted signal (10 MHz sine wave)
transmitted_signal = np.zeros_like(t)
transmitted_signal[:pulse_samples] = 10 * np.sin(2 * np.pi * f * t[:pulse_samples]) # Amplitude 10V
# Calculate round trip time for array 1
round_trip_time_1 = 2 * distance_1 / speed_of_sound # Total time for the sound to travel to the object and back
# Simulate received signal for array 1 with delay
received_signal_1 = np.zeros_like(t)
received_delay_index_1 = int(round_trip_time_1 / T)
received_signal_1[received_delay_index_1:received_delay_index_1 + pulse_samples] = transmitted_signal[:pulse_samples]
# Add random noise to the received signal for array 1
noise_1 = np.random.normal(0, 1, received_signal_1.shape) # Gaussian noise
received_signal_1 += noise_1
# Add DC offset of 100V
received_signal_1 += 100
# Calculate round trip time for array 2
round_trip_time_2 = 2 * distance_2 / speed_of_sound # Total time for the sound to travel to the object and back
# Simulate received signal for array 2 with delay
received_signal_2 = np.zeros_like(t)
received_delay_index_2 = int(round_trip_time_2 / T)
received_signal_2[received_delay_index_2:received_delay_index_2 + pulse_samples] = transmitted_signal[:pulse_samples]
# Add random noise to the received signal for array 2
noise_2 = np.random.normal(0, 1, received_signal_2.shape) # Gaussian noise
received_signal_2 += noise_2
# Add DC offset of 100V
received_signal_2 += 100
# Design Chebyshev I high-pass filter to remove DC component
nyquist = 0.5 * sampling_rate
cutoff_high = 1e5 # High-pass cutoff frequency (100 kHz)
order = 1
b, a = cheby1(order, 0.5, cutoff_high / nyquist, btype='highpass')
filtered_signal_high_1 = lfilter(b, a, received_signal_1)
filtered_signal_high_2 = lfilter(b, a, received_signal_2)
# Design Chebyshev I low-pass filter to remove high-frequency noise
cutoff_low = 15e6 # Low-pass cutoff frequency (15 MHz)
b_low, a_low = cheby1(order, 0.5, cutoff_low / nyquist, btype='low')
filtered_signal_low_1 = lfilter(b_low, a_low, filtered_signal_high_1)
filtered_signal_low_2 = lfilter(b_low, a_low, filtered_signal_high_2)
# Plot received and filtered signals for both arrays
plt.figure(figsize=(12, 10))
plt.subplot(4, 1, 1)
plt.title('Received Signal Array 1 with Noise and DC Offset')
plt.plot(t, received_signal_1, label='Received Signal Array 1', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 2)
plt.title('Received Signal Array 2 with Noise and DC Offset')
plt.plot(t, received_signal_2, label='Received Signal Array 2', color='blue')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 3)
plt.title('Filtered Signal Array 1 (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low_1, label='Filtered Signal Array 1', color='green')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 4)
plt.title('Filtered Signal Array 2 (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low_2, label='Filtered Signal Array 2', color='red')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.tight_layout()
plt.show()
# Cross-correlation for both arrays
correlation_1 = correlate(filtered_signal_low_1, transmitted_signal, mode='full')
correlation_2 = correlate(filtered_signal_low_2, transmitted_signal, mode='full')
lag = np.arange(-len(transmitted_signal) + 1, len(filtered_signal_low_1))
# Find the delays
max_corr_index_1 = np.argmax(correlation_1)
max_delay_1 = lag[max_corr_index_1] * T
max_corr_index_2 = np.argmax(correlation_2)
max_delay_2 = lag[max_corr_index_2] * T
# Calculate distances based on delays
calculated_distance_1 = max_delay_1 * speed_of_sound / 2 # Distance to object from array 1
calculated_distance_2 = max_delay_2 * speed_of_sound / 2 # Distance to object from array 2
print(f'Maximum correlation delay for Array 1: {max_delay_1:.6f} seconds, Calculated Distance: {calculated_distance_1:.2f} meters')
print(f'Maximum correlation delay for Array 2: {max_delay_2:.6f} seconds, Calculated Distance: {calculated_distance_2:.2f} meters')运行结果：
Maximum correlation delay for Array 1: 0.020000 seconds, Calculated Distance: 10.00 meters Maximum correlation delay for Array 2: 0.020100 seconds, Calculated Distance: 10.05 meters
六、【待优化】多阵元多目标探测(以2阵元2目标为例)
信号生成：生成线性调频（LFM）脉冲信号。
信号接收：模拟两个阵元接收两个目标的信号，计算每个阵元的延迟。
噪声和偏移：为接收到的信号添加高斯噪声和直流偏移。
滤波：使用切比雪夫高通和低通滤波器去除直流成分和高频噪声。
交叉相关：对过滤后的信号进行交叉相关，识别信号中的峰值。
距离计算：根据检测到的峰值计算目标的距离。
运行这个程序将会模拟两个阵元对两个目标的检测，并输出检测到的距离。调整参数可以模拟不同的场景。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import correlate, cheby1, lfilter, find_peaks
# Constants
f0 = 10e6 # Start frequency of the LFM pulse (10 MHz)
distance_vertical = 10 # Distance to the first object in meters
distance_horizontal = 50 # Distance to the second object in meters
speed_of_sound = 1000 # Speed of sound in m/s
sampling_rate = 100e6 # Sampling rate (100 MHz)
# Time settings
pulse_duration = 10 / f0 # Duration of the pulse in seconds
pulse_samples = int(sampling_rate * pulse_duration) # Number of samples for the pulse
# Generate time array
T = 1 / sampling_rate
t = np.arange(0, 2 * pulse_duration + 0.1, T) # Time array for sending and receiving
# Generate LFM pulse
k = (1e10) # Chirp rate
lfm_pulse = np.zeros_like(t)
lfm_pulse[:pulse_samples] = np.sin(2 * np.pi * (f0 * t[:pulse_samples] + 0.5 * k * t[:pulse_samples]**2)) # LFM pulse
# Simulate received signals with delays for two arrays
received_signal_1 = np.zeros_like(t)
received_signal_2 = np.zeros_like(t)
# Calculate delays for both objects for both arrays
received_delay_index1_array1 = int((2 * distance_vertical) / speed_of_sound / T)
received_delay_index2_array1 = int((2 * distance_horizontal) / speed_of_sound / T)
received_delay_index1_array2 = int((2 * distance_vertical) / speed_of_sound / T)
received_delay_index2_array2 = int((2 * distance_horizontal) / speed_of_sound / T)
# Simulate received signals
received_signal_1[received_delay_index1_array1:received_delay_index1_array1 + pulse_samples] = lfm_pulse[:pulse_samples]
received_signal_1[received_delay_index2_array1:received_delay_index2_array1 + pulse_samples] += lfm_pulse[:pulse_samples]
received_signal_2[received_delay_index1_array2:received_delay_index1_array2 + pulse_samples] = lfm_pulse[:pulse_samples]
received_signal_2[received_delay_index2_array2:received_delay_index2_array2 + pulse_samples] += lfm_pulse[:pulse_samples]
# Add random noise to the received signals (lower standard deviation for realistic simulation)
noise_std_dev = 0.1 # Standard deviation of noise
noise_1 = np.random.normal(0, noise_std_dev, received_signal_1.shape) # Gaussian noise
noise_2 = np.random.normal(0, noise_std_dev, received_signal_2.shape) # Gaussian noise
received_signal_1 += noise_1
received_signal_2 += noise_2
# Add DC offset of 100V
received_signal_1 += 100
received_signal_2 += 100
# Design Chebyshev I high-pass filter to remove DC component
nyquist = 0.5 * sampling_rate
cutoff_high = 1e5 # High-pass cutoff frequency (100 kHz)
order = 1
b, a = cheby1(order, 0.5, cutoff_high / nyquist, btype='highpass')
filtered_signal_high_1 = lfilter(b, a, received_signal_1)
filtered_signal_high_2 = lfilter(b, a, received_signal_2)
# Design Chebyshev I low-pass filter to remove high-frequency noise
cutoff_low = 15e6 # Low-pass cutoff frequency (15 MHz)
b_low, a_low = cheby1(order, 0.5, cutoff_low / nyquist, btype='low')
filtered_signal_low_1 = lfilter(b_low, a_low, filtered_signal_high_1)
filtered_signal_low_2 = lfilter(b_low, a_low, filtered_signal_high_2)
# Plot received and filtered signals
plt.figure(figsize=(12, 10))
plt.subplot(4, 1, 1)
plt.title('Received Signal Array 1 with Noise and DC Offset')
plt.plot(t, received_signal_1, label='Received Signal Array 1', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 2)
plt.title('Received Signal Array 2 with Noise and DC Offset')
plt.plot(t, received_signal_2, label='Received Signal Array 2', color='blue')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 3)
plt.title('Filtered Signal Array 1 (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low_1, label='Filtered Signal Array 1', color='green')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(4, 1, 4)
plt.title('Filtered Signal Array 2 (High-Pass and Low-Pass)')
plt.plot(t, filtered_signal_low_2, label='Filtered Signal Array 2', color='red')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.tight_layout()
plt.show()
# Cross-correlation for both arrays
correlation_1 = correlate(filtered_signal_low_1, lfm_pulse, mode='full')
correlation_2 = correlate(filtered_signal_low_2, lfm_pulse, mode='full')
lag = np.arange(-len(lfm_pulse) + 1, len(filtered_signal_low_1))
# Identify peaks in cross-correlation for both arrays
peaks_1, _ = find_peaks(correlation_1)
peaks_2, _ = find_peaks(correlation_2)
# Threshold for peak detection
threshold_1 = 0.8 * np.max(correlation_1)
threshold_2 = 0.8 * np.max(correlation_2)
valid_peaks_1 = peaks_1[correlation_1[peaks_1] > threshold_1]
valid_peaks_2 = peaks_2[correlation_2[peaks_2] > threshold_2]
# Calculate distances for detected objects
detected_distances_1 = []
detected_distances_2 = []
tolerance = 0.01 # Tolerance for detecting unique distances
last_distance_1 = None
last_distance_2 = None
for peak in valid_peaks_1:
distance_detected = abs((peak - len(lfm_pulse) + 1) * T * speed_of_sound / 2)
if last_distance_1 is None or abs(distance_detected - last_distance_1) > tolerance:
detected_distances_1.append(distance_detected)
last_distance_1 = distance_detected
for peak in valid_peaks_2:
distance_detected = abs((peak - len(lfm_pulse) + 1) * T * speed_of_sound / 2)
if last_distance_2 is None or abs(distance_detected - last_distance_2) > tolerance:
detected_distances_2.append(distance_detected)
last_distance_2 = distance_detected
print(f'Detected distances for Array 1: {detected_distances_1}')
print(f'Detected distances for Array 2: {detected_distances_2}')
计算结果：
Detected distances for Array 1: [9.9999, 49.9999] Detected distances for Array 2: [9.99995, 49.99995]
PS：其实对于干扰不是特别强烈的互相关计算来说，不进行低通滤波也是可以的，尤其是对于FPGA来说，这点很重要。
七、【待优化】单阵元单目标的FPGA实现(Xilinx HLS)
#include <hls_stream.h>
#include <ap_int.h>
#include <cmath>
#include <iostream>
#include <complex>
#include <cstdlib>
#define N 1024 // FFT Size
#define PULSE_SAMPLES 1000 // Number of samples for the pulse
#define SAMPLING_RATE 100e6
#define F 10e6
#define DISTANCE 10
#define SPEED_OF_SOUND 1000
typedef std::complex<float> complex_float;
// Function to generate transmitted signal
void generate_signal(float transmitted_signal[PULSE_SAMPLES]) {
for (int i = 0; i < PULSE_SAMPLES; i++) {
transmitted_signal[i] = 10 * sin(2 * M_PI * F * (i / SAMPLING_RATE));
}
}
// Function to simulate received signal with delay and noise
void simulate_received_signal(float transmitted_signal[PULSE_SAMPLES], float received_signal[N]) {
float round_trip_time = 2 * DISTANCE / SPEED_OF_SOUND;
int received_delay_index = static_cast<int>(round_trip_time * SAMPLING_RATE);
// Initialize received signal
for (int i = 0; i < N; i++) {
received_signal[i] = 0;
}
// Simulate received signal with delay
for (int i = 0; i < PULSE_SAMPLES; i++) {
if (i + received_delay_index < N) {
received_signal[i + received_delay_index] = transmitted_signal[i] + (rand() % 100) / 100.0; // Add noise
}
}
}
// Function to apply a simple high-pass filter
void high_pass_filter(float input[N], float output[N], float cutoff) {
float previous = 0;
float alpha = 1.0 / (1.0 + (SAMPLING_RATE / (2 * M_PI * cutoff))); // Simple RC filter
for (int i = 0; i < N; i++) {
output[i] = alpha * (input[i] - previous);
previous = input[i];
}
}
// Function to perform FFT (Cooley-Tukey)
void fft(complex_float input[N], complex_float output[N]) {
// Bit-reversal permutation
int n = N;
int j = 0;
for (int i = 0; i < n; i++) {
if (i < j) {
std::swap(input[i], input[j]);
}
int m = n / 2;
while (m >= 1 && j >= m) {
j -= m;
m /= 2;
}
j += m;
}
// FFT computation
for (int s = 1; s <= log2(n); s++) {
int m = 1 << s;
complex_float wm = std::exp(complex_float(0, -2 * M_PI / m));
for (int k = 0; k < n; k += m) {
complex_float w(1);
for (int j = 0; j < m / 2; j++) {
complex_float t = w * input[k + j + m / 2];
complex_float u = input[k + j];
input[k + j] = u + t;
input[k + j + m / 2] = u - t;
w *= wm;
}
}
}
for (int i = 0; i < n; i++) {
output[i] = input[i];
}
}
// Function to perform cross-correlation
void cross_correlation(float signal1[N], float signal2[PULSE_SAMPLES], float correlation[N]) {
for (int lag = -PULSE_SAMPLES + 1; lag < N; lag++) {
correlation[lag + PULSE_SAMPLES - 1] = 0;
for (int i = 0; i < PULSE_SAMPLES; i++) {
if (i + lag >= 0 && i + lag < N) {
correlation[lag + PULSE_SAMPLES - 1] += signal1[i + lag] * signal2[i];
}
}
}
}
// Top-level function
extern "C" void signal_processing(float transmitted_signal[PULSE_SAMPLES], float received_signal[N], float filtered_signal[N], float correlation[N]) {
// Step 1: Generate transmitted signal
generate_signal(transmitted_signal);
// Step 2: Simulate received signal
simulate_received_signal(transmitted_signal, received_signal);
// Step 3: High-pass filter
high_pass_filter(received_signal, filtered_signal, 1e5);
// Step 4: FFT
complex_float fft_input[N];
complex_float fft_output[N];
for (int i = 0; i < N; i++) {
fft_input[i] = filtered_signal[i];
}
fft(fft_input, fft_output);
// Step 5: Cross-correlation
cross_correlation(filtered_signal, transmitted_signal, correlation);
}
