超聲流量傳感器信號(hào)采樣率的影響

我們?cè)谥霸凇读髁?a href="http://wenjunhu.com/v/tag/117/" target="_blank">傳感器(2)—超聲流量傳感器，相位差和相關(guān)性原理》中曾留了一個(gè)問(wèn)題：如果采樣頻率不一樣，模擬的結(jié)果會(huì)怎么樣？不知道大家有沒(méi)有時(shí)間去運(yùn)行測(cè)試一下，反正小編嘗試了。測(cè)試的結(jié)果（但不是最終結(jié)論，請(qǐng)根據(jù)實(shí)際應(yīng)用進(jìn)行調(diào)整）下文會(huì)提到。

在小編上次關(guān)于超聲流量的相關(guān)信號(hào)處理的模擬中，用到了插值的方式，不過(guò)由于沒(méi)有經(jīng)過(guò)太多的測(cè)試，在后續(xù)測(cè)試過(guò)程中，發(fā)現(xiàn)當(dāng)聲速為4500m/s的模擬流體的流速不超過(guò)15m/s時(shí)，都沒(méi)有什么問(wèn)題;但是當(dāng)模擬的流速超過(guò)15m/s之后，兩組模擬信號(hào)的相位差超過(guò)了180°時(shí)，獲取相關(guān)性序列中對(duì)應(yīng)的時(shí)間值時(shí)，發(fā)現(xiàn)取值結(jié)果為負(fù)數(shù)。

相位差180°，這本身沒(méi)有問(wèn)題，但是由于數(shù)據(jù)的相關(guān)性處理函數(shù)是個(gè)偶函數(shù)，對(duì)稱(chēng)于中點(diǎn)，所以在數(shù)據(jù)相關(guān)性處理中，當(dāng)兩路信號(hào)的相位差超過(guò)180°前后時(shí)，因?yàn)閮陕窚y(cè)試信號(hào)是正余弦形式，相當(dāng)于其中一路信號(hào)序列在相對(duì)于另外一路信號(hào)序列移動(dòng)的時(shí)候，在某個(gè)相對(duì)處的相關(guān)計(jì)算結(jié)果，可以找到另外一個(gè)錯(cuò)位的對(duì)方有相同的相關(guān)處理結(jié)果，然而在查找標(biāo)識(shí)相位偏差所對(duì)應(yīng)的最大值時(shí)間序列值時(shí)，測(cè)試代碼優(yōu)先輸出了最左側(cè)編號(hào)是負(fù)數(shù)的相位序列，實(shí)際應(yīng)該輸出相關(guān)性序列中點(diǎn)右側(cè)的最大值所對(duì)應(yīng)的時(shí)間序列值。結(jié)果可想而知，計(jì)算的流速變成了負(fù)數(shù)。修改后的代碼中，計(jì)算流速時(shí)，使用的相關(guān)性序列是整個(gè)相關(guān)結(jié)果序列的右邊部分，就解決了輸出流速突然變成負(fù)數(shù)的問(wèn)題。

另外，實(shí)際的超聲波接收信號(hào)并不會(huì)是從頭至尾都是理想的正余弦，因此，我們看到的模擬信號(hào)似乎也是太過(guò)于理想。這里我們稍微加了點(diǎn)料，讓穩(wěn)定的信號(hào)后面加上了快速衰減。

在進(jìn)行模擬信號(hào)相關(guān)性處理的時(shí)候，之前我們使用的是直接死磕方式的計(jì)算，這里調(diào)整為快速傅里葉計(jì)算(FFT)。

所以，這里我們將模擬代碼作了以下調(diào)整，順便作了在不同采樣頻率下的模擬測(cè)量值比較：

在原來(lái)的模擬超聲波接收信號(hào)的生成數(shù)據(jù)中，加入最后衰減的部分;

在生成的相關(guān)性序列中，只取中間點(diǎn)右側(cè)的一半作為分析處理部分;

圖像輸出時(shí)，仍然輸出全部的相關(guān)性序列;

在生成相關(guān)性序列時(shí)，需用了FFT的模式(原先的選用的直算的方式);

降低“采樣頻率”到5MHz，然后采樣插值的方式進(jìn)行相關(guān)性處理。

上圖模擬的是16m/s下，5MHz采樣率，得到的輸出值：

T1_0: 4.5617943222309144e-07

Length of t lags: 2000

Ts Delayed cycles: 23 ,

Delayed time: 4.6e-07

degree(相位差，角度): 165.6

Sample time interval(采樣間隔時(shí)間): 2.0e-08

Estimated velocity: 16.134001423268288(m/s)

不同采樣率，不同模擬流速下的模擬輸出比較（都有插值）

從模擬結(jié)果來(lái)看，較高的采樣率確實(shí)可以獲得較高的測(cè)量精度，而低的采樣率，雖然適用了插值提高了測(cè)量精度，但是插值運(yùn)算畢竟是基于我們?cè)O(shè)定了波形的形態(tài)進(jìn)行的。實(shí)際應(yīng)用，還是要根據(jù)實(shí)際需求和情況配置相關(guān)的軟硬件功能。

另外需要說(shuō)明的是，在整個(gè)信號(hào)相關(guān)性處理的模擬過(guò)程中，這里并沒(méi)有將相關(guān)性結(jié)果序列歸一化后輸出，而是使用直接結(jié)果的方式。

小編還要留一個(gè)問(wèn)題：如果插值處理中標(biāo)識(shí)插值點(diǎn)數(shù)的參數(shù)有變化的時(shí)候，模擬結(jié)果會(huì)出現(xiàn)什么變化?歡迎大家測(cè)試并留個(gè)回復(fù)。

修改之后的模擬代碼

其中：模擬代碼中的噪聲部分是被注釋掉的

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate
import math
from scipy.interpolate import interp1d


def cor_demo():
    try:
        # Parameters
        c = 4500    # Sound speed in the fluid in m/s
        d = 0.5     # Pipe diameter in meters
        theta = math.pi / 6     # Angle of the emitted ultrasound wave, in radians. 30 degrees here.
        L = d / math.cos(theta) # The sound path length in meters parallel to the flow direction
        v = 16                  # Fluid velocity in m/s


        f = 1e6     # Frequency of the signal in 1MHz
        T = 1 / f   # Period of the signal in s


        Fs = 5e6        # 采樣頻率為5MHz
        Ts = 1/Fs       # 采樣間隔
        
        # Calculate time delay caused by flow velocity
        time_up = L / (c - v * math.sin(theta))     # Time of flight upstream
        time_down = L / (c + v * math.sin(theta))   # Time of flight downstream
        
        print("T1_0:",time_up-time_down)
        t_delay = 2 * L * v * np.sin(theta) / (c**2 - v**2 * np.sin(theta)**2)
        
        # Time array
        # 返回一個(gè)從0到period*T范圍內(nèi)的有num_samples個(gè)元素的一維數(shù)組，數(shù)組中的數(shù)值是等間距分布的
        period = 40
        num_samples = period*T/Ts
        t = np.linspace(0, period*T, int(num_samples))
        
        t_max = t.max()
        t_start_decay = t_max * 5/6


        # 定義信號(hào)衰減準(zhǔn)則 - 衰減系數(shù)一般由物質(zhì)的特性決定
        attenuation_coefficient = 10e5
        
        # Signals
        stationary_s0 = np.sin(2 * np.pi * f * t[t=t_start_decay])* np.exp(-attenuation_coefficient * (t[t >= t_start_decay] - t_start_decay))
        decay_s1 = np.sin(2 * np.pi * f * (t[t>=t_start_decay]-t_delay))* np.exp(-attenuation_coefficient * (t[t >= t_start_decay] - t_start_decay))
        
        #s0 = np.sin(2 * np.pi * f * t)* np.exp(-attenuation_coefficient * t)
        #s1 = np.sin(2 * np.pi * f * (t - t_delay))* np.exp(-attenuation_coefficient * t)  #*2.0 # Signal 1
        
        s0 = np.concatenate([stationary_s0, decay_so])
        s1 = np.concatenate([stationary_s1, decay_s1])
        """
        noise0 = np.random.normal(0, 0.5, s0.shape)
        noise1 = np.random.normal(0, 0.5, s1.shape)
        
        s0 = s0 + noise0
        s1 = s1 + noise1
        """
        # Define interpolation factor
        interp_factor = 10


        # New time vector after interpolation
        t_new = np.linspace(t.min(), t.max(), t.size * interp_factor)
        
        # Create a function based on the original signals, which can be used to generate the interpolated signals
        interp_func_s0 = interp1d(t, s0, kind='cubic')
        interp_func_s1 = interp1d(t, s1, kind='cubic')


        # Generate the interpolated signals
        s0_new = interp_func_s0(t_new)
        s1_new = interp_func_s1(t_new)


        # 計(jì)算Cross-correlation，并找到最大值對(duì)應(yīng)的位置
        #correlation = correlate(s1, s0, method='direct', mode='full')  # old
        correlation = correlate(s1_new, s0_new, method='fft', mode='full')
        
        # 找出該序列的長(zhǎng)度 len_corr 和其中間點(diǎn) mid_index
        len_corr = len(correlation)
        mid_index = len_corr // 2
        
        #取序列中間值右邊的序列（包含中間值）
        corr_half = correlation[mid_index:]


        #lags = np.arange(-len(s1_new) + 1, len(s1_new)) # Lags array
        lags = np.arange(0, len(s1_new)) # Lags array
        
        print("Length of t lags:", len(lags))
        # Calculate flow speed using the estimated time delay
        # Find the peak of the cross-correlation corresponds to the time delay
        delay = lags[np.argmax(corr_half)]            # 相位差所對(duì)應(yīng)的信號(hào)序列值
        
        # 采樣時(shí)間
        sample_time = (period*T) / len(t_new)
        
        time_delay = delay * sample_time                # 兩個(gè)信號(hào)間的延遲時(shí)間
        
        phase_shift = (time_delay / T) * 2 * np.pi      # 由時(shí)間延遲換算成的兩個(gè)信號(hào)序列的相位差
        phase_shift_deg = phase_shift * (180 / np.pi)   # 由相位差換成的兩個(gè)信號(hào)序列的角度差
        
        print("Ts Delayed cycles:", delay, ", Delayed time: ", time_delay, "degree:", phase_shift_deg) # 將延遲轉(zhuǎn)換為相位差
        print("Sample time interval:",sample_time)
        # 計(jì)算所得的流速
        v_estimated = (math.sqrt((L**2)+(time_delay**2)*(c**2))-L)/(time_delay * math.sin(theta))
        print('Estimated velocity: ', v_estimated)


        # 輸出圖
        fig, (ax_origin, ax_interpo, ax_corr) = plt.subplots(3, 1, figsize=(12, 8))


        # 原始信號(hào)圖
        ax_origin.plot(t, s0, label='Signal 1')
        ax_origin.plot(t, s1, label='Signal 2')
        ax_origin.set_title('Original signals')
        ax_origin.legend()


        # 插值后的信號(hào)圖
        ax_interpo.plot(t_new, s0_new, label='Signal 1')
        ax_interpo.plot(t_new, s1_new, label='Signal 2')
        ax_interpo.set_title('interpolated signals')
        ax_interpo.legend()
        # 互相關(guān)圖
        ax_corr.plot(correlation)
        ax_corr.axvline(x = len(correlation)//2 + delay, color = 'r', linestyle = '--', label = "Max correlation at delay")
        ax_corr.set_title('Cross-correlation between signal 1 and signal 2')
        ax_corr.legend()


        plt.tight_layout()
        plt.show()
        return
    except Exception as e:
        print("Error:",e)


if __name__=='__main__':
    cor_demo()