ultrasonic_m_echo_sin.m

在matlab上对复倒谱回波检测系统的仿真

clear all;close all;
%=============================
    fs=10E7;      % 采样率100MHz
    fTR=1E7;     % 发射波/回波中心频率10MHz
    WTR=fTR*2*pi; % 发射波/回波中心角频率,rad/S
    NT=1024;   
    NR=NT;   % 接收回波：50uS * 10MHz --> 500点
    Cus=6E3;
    Ts=1/fs;
    
    NoiseT=randn(1,NT);   % 产生随机噪声信号 
    tn=0:(NT-1);
    phi0=0.35*pi;
    vt=sin(WTR*tn*Ts+phi0) +0.001*NoiseT;  % 连续发射调频波
    kc=2/Cus/Ts;
    i=1;
    R1=8E-3;% 当调制波频率为5KHz/10KHz时，可达Rmax=50mm/35mm
    R1ind=round(R1*kc);
    vr1=[zeros(1,R1ind)  0.01*vt(1:(NT-R1ind))]; 
    R2=11E-3;% 当调制波频率为5KHz/10KHz时，可达Rmax=50mm/35mm
    R2ind=round(R2*kc);
    vr2=[zeros(1,R2ind)  0.008*vt(1:(NT-R2ind))]; 
    NoiseR=randn(1,NR);   % 产生随机噪声信号
    for R0=(3:0.1:12)*1E-3;% 当调制波频率为5KHz/10KHz时，可达Rmax=50mm/35mm
        R0ind=round(R0*kc)
        vr0=[zeros(1,R0ind)  0.2*vt(1:(NT-R0ind))]; 
        vr=vr0+vr1+vr2+0.001*NoiseR; % 叠加噪声, SNR~ = [0.04/0.01 0.025/0.01] = [4 2.5] =[12dB 8dB]
        v=vt+vr;
        vv(i,1:NT)=v;
        Vccep(i,1:NT)=abs(cceps(v));
        i=i+1;
    end

    Strind=80;
    Stpind=500;
    for i=1:length(R0)
        [Vmax(i) Vind(i)]=max(Vccep(i,Strind:Stpind)); % 实际回波位置范围
    end
    Rind=Vind+Strind-2
% % % % % % % % % % % %     VF=fft(v); 
% % % % % % % % % % % %     NVF=length(VF);  % 要求是偶数
% % % % % % % % % % % %     AngVF=phase(VF);
% % % % % % % % % % % %     LogVF=log(abs(VF));
% % % % % % % % % % % %     rpi=AngVF(NVF/2+1);`
% % % % % % % % % % % %     V=LogVF+1j*(AngVF-2*rpi/NVF);
% % % % % % % % % % % %     if LogVF(1)<0
% % % % % % % % % % % %         signV=-1;
% % % % % % % % % % % %     else
% % % % % % % % % % % %         signV=1;
% % % % % % % % % % % %     end
% % % % % % % % % % % %     
% % % % % % % % % % % %     vev=ifft(signV*V);
    
% Differentiator design - Hamming window    
% % % % % %     M = 11; alpha = (M-1)/2;
% % % % % %     n = 0:M-1;
% % % % % %     hd = (cos(pi*(n-alpha)))./(n-alpha+eps);
% % % % % %     hd(alpha+1)=0;
% % % % % %     w_ham = (hamming(M))';
% % % % % %     hdif = hd .* w_ham;
% % % % % %     [Hr,w,P,L] = hr_type3(hdif);
% % % % % %     y=abs(filter(hdif, 1, vev));
    

    N=21;
    f=[0:0.05:1-0.05];
    m=10*f;
    [hdif] = REMEZ(N, f, m, 'differentiator');
    y=abs(filter(hdif, 1, Vccep));
% % % % % [h] = REMEZ(N, f, m, weight, ftype)
% % % % % N为阶数，h长度为：a = N+1
% % % % % f为边界频率，从小到大的次序
% % % % % m为与边界频率f对应的幅度，m，f同长，同偶
% % % % % weight为与频带对应的权值，缺省权值为1
% % % % % ftype为滤波器类型：
% % % % % 	‘hilbert‘，  ’differentiator’，
% % % % % 	缺省为选频滤波器
    
    
% % %     VT=hilbert(vt);
% % %     VR=hilbert(vr);
% % %     phiT=phase(VT);
% % %     phiR=phase(VR);
% % %     deltaphi=phiR-phiT;
    
    figure(1);
    plot(Vccep(15,1:NT));
    figure(2);
    mesh(Vccep);
    grid;
    figure(3);
    mesh(y);
    grid;


% % %     hold on;
% % %     h=plot(phiR,'r--');
% % %     set(h,'linewidth',2);
% % %     grid;
% % %     hold off;
% % %     subplot(212)
% % %     h=plot(deltaphi,'r');
% % %     set(h,'linewidth',2);
% % %     grid;
% % %     deltaphiav=mean(deltaphi(30:length(deltaphi)-30));
% % %     set(gca,'XTickMode','manual','YTick',deltaphiav)
% % % 
% % %     Rind=0.5*deltaphiav*Cus/WTun
% % %     Range=1E3*Rind