FPGA-Matlab信号处理笔记

1 ILA读取IQ信号

1.1 解析

ILA抓取IQ信号且valid不连续，或实现抽取时数据不连续，此时用valid判断再读出valid；

将IQ转为实信号，使用complex函数。

将数据读出到工作区时，可直接使用内APP “Signal Analyzer”分析时频。

1.2 代码实现

clc;

clear all;

close all;

im_i = readmatrix('iladata_ad.csv', 'outputtype', 'string'); % 读取string矩阵

re_q = readmatrix('iladata_ad.csv', 'outputtype', 'string'); % 读取string矩阵

% indexs=(double(im_i(:,8))==-1);

im_ddc = im_i((double(im_i(:,10))==1),12); % 取第x列

re_ddc = re_q((double(re_q(:,10))==1),11); % 取第x列

ddc_data_i = str2double(im_ddc);

ddc_data_q = str2double(re_ddc);

data_real = complex(ddc_data_i,ddc_data_q);

fs = 250*1e6 ;

dec =1;

data_iq = ddc_data_i + 1j * ddc_data_q ;

Plot_Freq_Jason(data_iq.',fs/dec,1,'DDC_Plot')

2 插值抽取

2.1 解析说明

插4抽5实现采样率250M降至200M。

2.2 代码实现

clear all;

close all;

k = 4; % 测试插入倍数

fc4 = 200 *1e6; % 2.5 * 4 = 10M% 插值后的采样率

% 插值后的滤波器系数

Fs = 2; % Sampling Frequency

Fpass = 0.10; % Passband Frequency  %%注意带宽

Fstop = 0.20; % Stopband Frequency

Dpass = 0.0057501127785; % Passband Ripple

Dstop = 0.000000001; % Stopband Attenuation

dens = 20; % Density Factor

[N, Fo, Ao, W] = firpmord([Fpass, Fstop] / (Fs / 2), [1 0], [Dpass,Dstop]);

 fir_n = firpm(N, Fo, Ao, W, {dens});

% 多相滤波器系数 一分四

fir_len = length(fir_n);

fir_n_1 = fir_n(1:4:fir_len);

fir_n_2 = fir_n(2:4:fir_len);

fir_n_3 = fir_n(3:4:fir_len);

fir_n_4 = fir_n(4:4:fir_len);

s = ddc_data_i ;%生成/读取一组AD信号--抽取之前的250M信号

% 多相插值滤波(k = 4)

y_1 = filter(fir_n_1, 1, s) * k;

y_2 = filter(fir_n_2, 1, s) * k;

y_3 = filter(fir_n_3, 1, s) * k;

y_4 = filter(fir_n_4, 1, s) * k;

y_n = [y_1; y_2; y_3; y_4];

yy = y_n(:)';%插值完成

yy1 = yy(1:5:end);%抽5

Plot_Freq_Jason(yy1.',fc4,1,'DDC_Plot')

3 SFNR计算

无杂散动态范围

N = length(data_real); %读取原始AD数据

f = (0:N-1) * (Fs/N)    ; % 频率轴

X= fft(data_real); %信号fft

P = abs(X/N);%功率谱

%计算SFDR : 无杂散动态范围

[~,maxSigidx]  = max(P(1:N));

maxSigpower = P(maxSigidx); %最大信号功率

P(maxSigidx-100:maxSigidx+100) = 0 ; %屏蔽最大信号周围频率

[~,maxSpurIdx] = max(P(1:N/2+1)); %最大信号功率

maxSpurpower = P(maxSpurIdx); %最大杂散功率

SFDR = maxSigpower - maxSpurpower ;

SFDR_dB = 20*log10(SFDR);%转换成dB

fprintf('SFDR =  %ddB\n',SFDR_dB) ;

4 相位计算

4.1 说明

将两通道的数据过拟合成y = a + b * sin(2 * pi * d * t + c)

可得出幅度b、初相c、频率d、直流a；

Δc= (c1 - c2)/pi*180°/180/d；相位差单位us

求取相位差总体思路是：得到两通道的频谱数据，且信号频率一致，使用最小二乘法转换成正弦公式，再求取相位差数据。

4.2 代码实现

% 文件读取

ch1 = adc_ch1_data;%自行导入通道数据

ch2 = adc_ch2_data;%自行导入通道数据

fs = 1e6; %采样率

% 初步比较信号是不是一致的

fft_c1 = abs(fft(ch1)); %对两通道数据做FFT

fft_c2 = abs(fft(ch2)); %对两通道数据做FFT

[~, b1] = max(fft_c1(2:end)); %找最大值的位置

[~, b2] = max(fft_c2(2:end)); %找最大值的位置

if (b1 ~= b2)

    error("error: f1 != f2!");

end

% 估算信号上下限

ff = fs / length(ch1);

f_b = b1 * ff; %找到信号频率

f1 = max(0, f_b - ff); %细化起始频率

f2 = f_b + ff; %细化结束频率

% 最小二乘法 拟合y = a + b * sin(2 * pi * d * t + c);

[a1, b1, c1, d1, e1] = sin_t(ch1, fs, f1, f2);

[a2, b2, c2, d2, e2] = sin_t(ch2, fs, f1, f2);

fprintf("CH1 直流:%f,幅值:%f,初相:%f°,频率:%f,标准差:%f\n", a1, b1, c1 / pi * 180, d1, e1);

fprintf("CH2 直流:%f,幅值:%f,初相:%f°,频率:%f,标准差:%f\n", a2, b2, c2 / pi * 180, d2, e2);

% 显示原始数据和拟合数据

t = (0:(length(ch1) - 1)) / fs;

ch1_t = a1 + b1 * cos(2 * pi * d1 * t + c1);

ch2_t = a2 + b2 * cos(2 * pi * d2 * t + c2);

% 计算两通道的相位差

pha_ba_diff = (c1 - c2) / pi * 180;

if (pha_ba_diff < -180)

    pha_ba_diff = pha_ba_diff + 360;

elseif (pha_ba_diff > 180)

    pha_ba_diff = pha_ba_diff - 360;

end

delay_t_diff = pha_ba_diff / 180 / d1; %相位差计算结果，单位为us

if (abs(delay_t_diff) < 1e-9)

    delay_t_diff = delay_t_diff * 1e12;

    delay_t_unit = 'ps';

elseif (abs(delay_t_diff) < 1e-6)

    delay_t_diff = delay_t_diff * 1e9;

    delay_t_unit = 'ns';

elseif (abs(delay_t_diff) < 1e-3)

    delay_t_diff = delay_t_diff * 1e6;

    delay_t_unit = 'us';

elseif (abs(delay_t_diff) < 1)

    delay_t_diff = delay_t_diff * 1e3;

    delay_t_unit = 'ms';

else

    delay_t_unit = 's';

end

fprintf("CH1-CH2 时间差:%f %s\n", delay_t_diff, delay_t_unit);

fprintf("CH1-CH2 相位差:%f °\n", delay_t_diff/(10^12)/(1/1000000)*360);%根据计算的时间差和自身频率周期进行换算

%% functions

function [a,b,c,d,e]=sin_t(y, fs, fl, fh)

% y = a + b*cos(2*pi*d*t + c);

if(length(y(1,:)) ~= 1)

    y = y';

end

N = length(y);

t = ((0:(N-1))/fs)';

cfun = @(d) [ones(size(t)), sin(2*pi*d*t), cos(2*pi*d*t)]\y;

sumerr2 = @(d) sum((y - [ones(size(t)), sin(2*pi*d*t), cos(2*pi*d*t)] * cfun(d)) .^ 2);

dopt = fminbnd(sumerr2, fl, fh);

abb = cfun(dopt);

a = abb(1);

b = norm(abb([2 3]));

c = atan2(-abb(2), abb(3));%求arctan

d = dopt;

e = sqrt(sumerr2(dopt) / N);

end