位置: IT常识 - 正文

开源代码 | FMCW-MIMO雷达仿真MATLAB(开源代码网站github)

编辑:rootadmin
开源代码 | FMCW-MIMO雷达仿真MATLAB

推荐整理分享开源代码 | FMCW-MIMO雷达仿真MATLAB(开源代码网站github),希望有所帮助,仅作参考,欢迎阅读内容。

文章相关热门搜索词:开源代码在哪找,开源代码网站,开源代码公布意味着什么,开源代码与组件使用情况说明怎么写,开源代码商用是否违法,开源代码是什么意思,开源代码与组件使用情况说明怎么写,开源代码公布意味着什么,内容如对您有帮助,希望把文章链接给更多的朋友!

本文编辑:调皮哥的小助理

本程序来源:https://github.com/ekurtgl/FMCW-MIMO-Radar-Simulation,作者是阿拉巴马大学博士生艾库特格尔,研究方向主要是雷达信号处理人类活动识别以及雷达数据的机器学习应用,这份比较新的开源雷达仿真代码,值得大家学习。

下面主要分析代码的主要内容,方便大家解读。

程序目录如下:

图片

FMCW_simulation.m是创建点目标并估计其范围、速度和角度信息的主脚本,首先研究这个脚本程序。

一、雷达参数

雷达参数的设置,属于是老生常谈了,之前的文章已经谈到很多了,不再详细重复论述。不过在本程序中,需要注意PRI默认为等于Chirp持续时长,实际上还要考虑空闲时间。这里的带宽指的是有效带宽,而不是发射信号带宽。程序中一共设置了10帧、1T8R,这些参数都是可以修改的。

```c> %% Radar parametersc = physconst('LightSpeed'); %speed of lightBW = 150e6; %bandwidth 有效fc = 77e9; % carrier frequencynumADC = 256; % # of adc samplesnumChirps = 256; % # of chirps per framenumCPI = 10;T = 10e-6; % PRI,默认不存在空闲时间PRF = 1/T;Fs = numADC/T; % sampling frequencydt = 1/Fs; % sampling intervalslope = BW/T;lambda = c/fc;N = numChirps*numADC*numCPI; % total # of adc samplest = linspace(0,T*numChirps*numCPI,N); % time axis, one frame 等间隔时间/点数t_onePulse = 0:dt:dt*numADC-dt; %单chirp时间numTX = 1;numRX = 8; %等效后Vmax = lambda/(T*4); % Max Unamb velocity m/sDFmax = 1/2*PRF; % = Vmax/(c/fc/2); % Max Unamb Dopp FreqdR = c/(2*BW); % range resolRmax = Fs*c/(2*slope); % TI's MIMO Radar docRmax2 = c/2/PRF; % lecture 2.3dV = lambda/(2*numChirps*T); % velocity resol, lambda/(2*framePeriod)d_rx = lambda/2; % dist. between rxsd_tx = 4*d_rx; % dist. between txsN_Dopp = numChirps; % length of doppler FFTN_range = numADC; % length of range FFTN_azimuth = numTX*numRX;R = 0:dR:Rmax-dR; % range axisV = linspace(-Vmax, Vmax, numChirps); % Velocity axisang_ax = -90:90; % angle axis二、目标参数

这里的目标参数就采用极坐标系转直角坐标系,然后分别用X方向的速度和Y方向上的速度乘上时间,带入后续的计算。

> %% 目标参数r1_radial = 50; v1_radial = 10; % velocity 1tar1_angle = -10;r1_y = cosd(tar1_angle)*r1_radial;r1_x = sind(tar1_angle)*r1_radial;v1_y = cosd(tar1_angle)*v1_radial;v1_x = sind(tar1_angle)*v1_radial;r1 = [r1_x r1_y 0];r2_radial = 100;v2_radial = -15; % velocity 2tar2_angle = 10;r2_y = cosd(tar2_angle)*r2_radial;r2_x = sind(tar2_angle)*r2_radial;v2_y = cosd(tar2_angle)*v2_radial;v2_x = sind(tar2_angle)*v2_radial;r2 = [r2_x r2_y 0];%发射天线位置tx_loc = cell(1,numTX);for i = 1:numTX tx_loc{i} = [(i-1)*d_tx 0 0]; scatter3(tx_loc{i}(1),tx_loc{i}(2),tx_loc{i}(3),'b','filled') hold onend% 接收天线位置rx_loc = cell(1,numRX);for i = 1:numRX rx_loc{i} = [tx_loc{numTX}(1)+d_tx+(i-1)*d_rx 0 0]; scatter3(rx_loc{i}(1),rx_loc{i}(2),rx_loc{i}(3),'r','filled')endtar1_loc = zeros(length(t),3);tar1_loc(:,1) = r1(1) + v1_x*t;tar1_loc(:,2) = r1(2) + v1_y*t;tar2_loc = zeros(length(t),3);tar2_loc(:,1) = r2(1) + v2_x*t;tar2_loc(:,2) = r2(2) + v2_y*t;

同时,这里还绘制了发射天线和接收天线的位置图,蓝色为发射天线位置,红色为接收天线位置。这个位置信息后面会用于计算目标的时间延迟。

三、发射信号建模

利用收发天线的位置以及目标参数中雷达的位置信息,先求目标到雷达的2-范数(也就是空间中两点的直线距离),然后转化为目标的延迟时间τ,如此以来得到的信号模型精度更高!这样的方法,与往常采用的不同,可以说以前的方法有点粗暴,虽同为远场条件,但是没有考虑阵列的位置对信号的影响!

> %% TX siganldelays_tar1 = cell(numTX,numRX);delays_tar2 = cell(numTX,numRX);r1_at_t = cell(numTX,numRX);r2_at_t = cell(numTX,numRX);tar1_angles = cell(numTX,numRX);tar2_angles = cell(numTX,numRX);tar1_velocities = cell(numTX,numRX);tar2_velocities = cell(numTX,numRX);for i = 1:numTX for j = 1:numRX delays_tar1{i,j} = (vecnorm(tar1_loc-repmat(rx_loc{j},N,1),2,2)+vecnorm(tar1_loc-repmat(tx_loc{i},N,1),2,2))/c; delays_tar2{i,j} = (vecnorm(tar2_loc-repmat(rx_loc{j},N,1),2,2)+vecnorm(tar2_loc-repmat(tx_loc{i},N,1),2,2))/c; endend四、接收信号模型

这里,接收信号没有采用发射与接收混频的形式,而是相位直接做差,分别计算两个目标的中频信号相加,此法等效为混频。其中@(tx,fx)相当于传参函数,用于后面的等式计算。

中频信号除了目标的距离表征的频率之外,还有目标运动速度引起的多普勒频移,是两部分的叠加。

> %% Complex signalphase = @(tx,fx) 2*pi*(fx.*tx+slope/2*tx.2); % transmittedphase2 = @(tx,fx,r,v) 2*pi*(2*fx*r/c+tx.*(2*fx*v/c + 2*slope*r/c)); % downconverted% f_oneChirp = slope*t(1:sum(t<=T));% f_t = repmat(f_oneChirp,1,numChirps*numCPI)-(BW/2); % transmit freq% f_t = BW/2*sawtooth(t/T*2*pi); fr1 = 2*r1(2)*slope/c; fr2 = 2*r2(2)*slope/c;fd1 = 2*v1_radial*fc/c; % doppler freqfd2 = 2*v2_radial*fc/c;f_if1 = fr1 + fd1; % beat or IF freqf_if2 = fr2 + fd2;% mixed1 = cell(numTX,numRX);% mixed2 = cell(numTX,numRX);mixed = cell(numTX,numRX);for i = 1:numTX for j = 1:numRX disp(['Processing Channel: ' num2str(j) '/' num2str(numRX)]); for k = 1:numChirps*numCPI phase_t = phase(t_onePulse,fc); phase_1 = phase(t_onePulse-delays_tar1{i,j}(k*numADC),fc); % received phase_2 = phase(t_onePulse-delays_tar2{i,j}(k*numADC),fc); signal_t((k-1)*numADC+1:k*numADC) = exp(1j*phase_t); signal_1((k-1)*numADC+1:k*numADC) = exp(1j*(phase_t - phase_1)); signal_2((k-1)*numADC+1:k*numADC) = exp(1j*(phase_t - phase_2)); end mixed{i,j} = signal_1 + signal_2; endend

绘制出的局部结果如下,若要观察全部的信号,则需要修改X轴的范围限制。

五、2D-FFT

上面我一句雷达信号处理原理方面的内容都没有提及,主要是这些基础的内容,我默认各位都已经清楚了。如果不清楚,那就先去搞清楚再回头来看(调皮连续波公众号里有),不然可能会遇到理解困难。

首先生成一个三维的DataCube,然后每个CPI做一次2D-FFT就OK了,没啥难度,顶多就是数据组成或者格式,理解起来有些麻烦,不过这都不是核心难题。

> RDC = reshape(cat(3,mixed{:}),numADC,numChirps*numCPI,numRX*numTX); % radar data cubeRDMs = zeros(numADC,numChirps,numTX*numRX,numCPI);for i = 1:numCPI RD_frame = RDC(:,(i-1)*numChirps+1:i*numChirps,:); RDMs(:,:,:,i) = fftshift(fft2(RD_frame,N_range,N_Dopp),2);endfigureimagesc(V,R,20*log10(abs(RDMs(:,:,1,1))/max(max(abs(RDMs(:,:,1,1))))));colormap(jet(256))% set(gca,'YDir','normal')clim = get(gca,'clim');caxis([clim(1)/2 0])xlabel('Velocity (m/s)');ylabel('Range (m)');

处理结果如下:

化成三维的可以标注数据,方便查看解算的结果以及误差,可见距离和速度都基本吻合,存在些许误差。

六、CA-CFAR开源代码 | FMCW-MIMO雷达仿真MATLAB(开源代码网站github)

代码一读就懂,老生常谈的东西了,没啥可讲的。

> numGuard = 2; % # of guard cellsnumTrain = numGuard*2; % # of training cellsP_fa = 1e-5; % desired false alarm rate SNR_OFFSET = -5; % dBRDM_dB = 10*log10(abs(RDMs(:,:,1,1))/max(max(abs(RDMs(:,:,1,1)))));[RDM_mask, cfar_ranges, cfar_dopps, K] = ca_cfar(RDM_dB, numGuard, numTrain, P_fa, SNR_OFFSET);figureh=imagesc(V,R,RDM_mask);xlabel('Velocity (m/s)')ylabel('Range (m)')title('CA-CFAR')

处理结果:

七、角度估计(一)3D-FFT

常规操作,基本内容,没啥难度。

```c> rangeFFT = fft(RDC(:,1:numChirps,:),N_range);angleFFT = fftshift(fft(rangeFFT,length(ang_ax),3),3);range_az = squeeze(sum(angleFFT,2)); % range-azimuth mapfigurecolormap(jet)% imagesc(ang_ax,R,20*log10(abs(range_az)./max(abs(range_az(:))))); mesh(ang_ax,R,20*log10(abs(range_az)./max(abs(range_az(:))))); xlabel('Azimuth Angle')ylabel('Range (m)')title('FFT Range-Angle Map')set(gca,'clim', [-35, 0])doas = zeros(K,181); % direction of arrivalsfigurehold on; grid on;for i = 1:K doas(i,:) = fftshift(fft(rangeFFT(cfar_ranges(i),cfar_dopps(i),:),181)); plot(ang_ax,10*log10(abs(doas(i,:))))endxlabel('Azimuth Angle')ylabel('dB')

运行结果如下: 可见,角度并不是很准确,读者可以验证,目标角度越偏离发现,误差越大。这其实与雷达的角度分辨率有关,用下文的超分辨算法可以改善。

(二)MUSIC算法

有读者常问,说MUSIC在网上找到的代码都是仿真的,很少看到有对实际目标的数据进行处理的,这不就来了嘛。不过,说到底还是见的太少。

MUSIC算法原理默认都清楚,不清楚的自己查阅相关资料。

> d = 0.5;M = numCPI; % # of snapshotsfor k=1:length(ang_ax) a1(:,k)=exp(-1i*2*pi*(d*(0:numTX*numRX-1)'*sin(ang_ax(k).'*pi/180)));endfor i = 1:K Rxx = zeros(numTX*numRX,numTX*numRX); for m = 1:M A = squeeze(RDMs(cfar_ranges(i),cfar_dopps(i),:,m)); Rxx = Rxx + 1/M * (A*A'); end [Q,D] = eig(Rxx); % Q: eigenvectors (columns), D: eigenvalues [D, I] = sort(diag(D),'descend'); Q = Q(:,I); % Sort the eigenvectors to put signal eigenvectors first Qs = Q(:,1); % Get the signal eigenvectors Qn = Q(:,2:end); % Get the noise eigenvectors for k=1:length(ang_ax) music_spectrum(i,k)=(a1(:,k)'*a1(:,k))/(a1(:,k)'*(Qn*Qn')*a1(:,k)); endend

运行结果漂亮的很,不要管幅度差异,幅度不代表功率。

(三)点云生成

坐标转换即可,若想要点云数量更多,可以降低CFAR门限,放更多的目标点进来。

> [~, I] = max(music_spectrum(2,:));angle1 = ang_ax(I);[~, I] = max(music_spectrum(1,:));angle2 = ang_ax(I);coor1 = [cfar_ranges(2)*sind(angle1) cfar_ranges(2)*cosd(angle1) 0];coor2 = [cfar_ranges(1)*sind(angle2) cfar_ranges(1)*cosd(angle2) 0];figurehold on;title('3D Coordinates (Point Cloud) of the targets')scatter3(coor1(1),coor1(2),coor1(3),100,'m','filled','linewidth',9)scatter3(coor2(1),coor2(2),coor2(3),100,'b','filled','linewidth',9)xlabel('Range (m) X')ylabel('Range (m) Y')zlabel('Range (m) Z')

处理结果:

(四)MUSIC 距离-AOA谱

算法没变,绘图方式不同而已。

> rangeFFT = fft(RDC);for i = 1:N_range Rxx = zeros(numTX*numRX,numTX*numRX); for m = 1:M A = squeeze(sum(rangeFFT(i,(m-1)*numChirps+1:m*numChirps,:),2)); Rxx = Rxx + 1/M * (A*A'); end% Rxx = Rxx + sqrt(noise_pow/2)*(randn(size(Rxx))+1j*randn(size(Rxx))); [Q,D] = eig(Rxx); % Q: eigenvectors (columns), D: eigenvalues [D, I] = sort(diag(D),'descend'); Q = Q(:,I); % Sort the eigenvectors to put signal eigenvectors first Qs = Q(:,1); % Get the signal eigenvectors Qn = Q(:,2:end); % Get the noise eigenvectors for k=1:length(ang_ax) music_spectrum2(k)=(a1(:,k)'*a1(:,k))/(a1(:,k)'*(Qn*Qn')*a1(:,k)); end range_az_music(i,:) = music_spectrum2;endfigurecolormap(jet)imagesc(ang_ax,R,20*log10(abs(range_az_music)./max(abs(range_az_music(:))))); xlabel('Azimuth')ylabel('Range (m)')title('MUSIC Range-Angle Map')clim = get(gca,'clim');

(五)压缩感知

原理捎带复杂,但主要还是一个凸优化问题。

> figurehold on; grid on;title('Angle Estimation with Compressed Sensing')xlabel('Azimuth')ylabel('dB')for i = 1:K A = squeeze(RDMs(cfar_ranges(i),cfar_dopps(i),:,1)); cvx_begin variable s(numTheta) complex; %alphax(numTheta,1) phix(numTX*numRX,numTheta)... %cap_theta(numTX*numRX,numTheta) %B(numTX*numRX,numTheta)%psix(numTheta,numTheta) %A(numRX*numTX,1) % A is the initial measurement % cap_theta == phix * psix; % minimize(norm(alphax,1)) % pow_p(norm(A-cap_theta*alphax,2),2) <= 1; % norm(A-cap_theta*alphax,2) <= 1; % minimize(norm(A-cap_theta*alphax,1)) minimize(norm(s,1)) norm(A-B*s,2) <= 1; cvx_end cvx_status cvx_optval plot(ang_ax,10*log10(abs(s)))end

第一个脚本的内容就说完了,下面是第二个脚本。

FMCW_sim_v2.m是读取 Kinect v2 设备捕获的人体骨骼关节三维坐标,提取原始雷达数据立方体 (RDC) 并播放距离多普勒图、输出微多普勒频谱图的主脚本。

Kinect v2里人体骨骼结构如下:

距离多普勒图:

微多普勒频谱图:

> rBin = 1:256;nfft = 212;window = 256;noverlap = 200;shift = window - noverlap;sx = myspecgramnew(sum(RDC(rBin,:,:)),window,nfft,shift); % mti filter and IQ correctionsx2 = abs(flipud(fftshift(sx,1)));timeAxis = [1:numCPI]*frameDuration; % TimefreqAxis = linspace(-PRF/2,PRF/2,nfft); % Frequency Axisfig = figure('visible','on');colormap(jet(256));% set(gca,'units','normalized','outerposition',[0,0,1,1]);doppSignMTI = imagesc(timeAxis,[-PRF/2 PRF/2],20*log10(abs(sx2/max(sx2(:)))));% axis xy% set(gca,'FontSize',10)title('micro-Doppler Spectrogram');% title(fOut(end-22:end-4))xlabel('Time (sec)');ylabel('Frequency (Hz)');caxis([-45 0]) % 40set(gca, 'YDir','normal')set(gcf,'color','w');

对该开源仿真代码感兴趣的读者,可以访问链接:https://github.com/ekurtgl/FMCW-MIMO-Radar-Simulation。

修改后的代码放到云盘了,修改的不多,如有需要,后台回复“2022”获得下载链接。

【本期结束】

本文链接地址:https://www.jiuchutong.com/zhishi/299086.html 转载请保留说明!

上一篇:【Tensorflow深度学习】实现手写字体识别、预测实战(附源码和数据集 超详细)(tensorflow gan)

下一篇:安装nvm,并使用nvm安装nodejs及配置环境变量(nvme安装win10教程)

  • 手机qq文件怎么重命名(手机qq文件怎么发到qq邮箱)

    手机qq文件怎么重命名(手机qq文件怎么发到qq邮箱)

  • 小米10安全模式怎么退出来啊(小米10安全模式后软件消失)

    小米10安全模式怎么退出来啊(小米10安全模式后软件消失)

  • 华为手机怎么硬格式化(华为手机怎么硬格机)

    华为手机怎么硬格式化(华为手机怎么硬格机)

  • 微信图片查看原图不显示了(微信图片查看原图是什么意思)

    微信图片查看原图不显示了(微信图片查看原图是什么意思)

  • 企业微信外部群20人怎么扩充(企业微信外部群人数上限)

    企业微信外部群20人怎么扩充(企业微信外部群人数上限)

  • 华为p30pro如何插耳机(华为p30 pro怎么连接u盘)

    华为p30pro如何插耳机(华为p30 pro怎么连接u盘)

  • 微信昵称能查到信息吗(微信昵称能查到什么信息)

    微信昵称能查到信息吗(微信昵称能查到什么信息)

  • iphone充电显示有液体(iphone充电显示有液体,然后无反应)

    iphone充电显示有液体(iphone充电显示有液体,然后无反应)

  • ascii码是几位码(ascii码是几位的)

    ascii码是几位码(ascii码是几位的)

  • 苹果x手机led灯怎么设置(苹果x设置led灯闪烁)

    苹果x手机led灯怎么设置(苹果x设置led灯闪烁)

  • 苹果突然出现请求登录(苹果手机出现请联系itunes是什么意思)

    苹果突然出现请求登录(苹果手机出现请联系itunes是什么意思)

  • 三星手机不断重启怎么回事(三星手机不断重启无法开机)

    三星手机不断重启怎么回事(三星手机不断重启无法开机)

  • 微信删除了好友还可以发信息吗(微信删除了好友聊天记录能恢复吗)

    微信删除了好友还可以发信息吗(微信删除了好友聊天记录能恢复吗)

  • 华为手机屏幕怎么显示步数(华为手机屏幕怎么设置)

    华为手机屏幕怎么显示步数(华为手机屏幕怎么设置)

  • 中文word是什么软件(world是什么)

    中文word是什么软件(world是什么)

  • 拼多多账号可以在两个手机上同时使用吗(拼多多账号可以在几个手机上登录)

    拼多多账号可以在两个手机上同时使用吗(拼多多账号可以在几个手机上登录)

  • 查找id对方有提醒吗(查找id定位对方)

    查找id对方有提醒吗(查找id定位对方)

  • 支付宝解除限制要多久(支付宝解除限制后多久可以转账)

    支付宝解除限制要多久(支付宝解除限制后多久可以转账)

  • 小米手机储存空间其他文件是什么(小米手机储存空间在哪里)

    小米手机储存空间其他文件是什么(小米手机储存空间在哪里)

  • 红外测温仪ems什么意思(红外测温仪ems什么牌子好)

    红外测温仪ems什么意思(红外测温仪ems什么牌子好)

  • 荣耀v30没有耳机孔怎么听歌(荣耀v30插耳机没反应)

    荣耀v30没有耳机孔怎么听歌(荣耀v30插耳机没反应)

  • 手机wps表格怎么画表(手机wps表格怎么调整表格大小)

    手机wps表格怎么画表(手机wps表格怎么调整表格大小)

  • 为什么我的华为p30pro充电慢(为什么我的华为mate40pro升级不了鸿蒙3.0)

    为什么我的华为p30pro充电慢(为什么我的华为mate40pro升级不了鸿蒙3.0)

  • 苹果耳机a1602是国行么(iphone耳机a1602)

    苹果耳机a1602是国行么(iphone耳机a1602)

  • logo设计说明怎么写(logo设计说明怎么写100字)

    logo设计说明怎么写(logo设计说明怎么写100字)

  • 安卓手机字体怎么改(安卓手机字体怎么改简体)

    安卓手机字体怎么改(安卓手机字体怎么改简体)

  • 个体户转到个人要多少税
  • 土地使用税返还是否征税
  • 什么是库存现金的盘亏
  • 公司组织旅游的费用要交个税
  • 红字发票开错了已上传如何作废
  • 哪些商业保险可以扣除个人所得税
  • 企业所得税申报表在哪里查询
  • 转让金融商品的会计分录
  • 旅行社差额征收怎么做账
  • 红字冲销增值税专用发票怎么写
  • 专利权转让的净收益计入
  • 营改增服务业税率
  • 企业如何规避印刷风险
  • 上报汇总和抄报是一个意思吗
  • 个人所得税特殊计税方法
  • 专项维修基金和契税有什么区别
  • 管家婆已过账销售单如何删除
  • 分支机构企业所得税是否必须跟总公司分摊吗
  • 建筑行业收到劳务发票入工程施工科目
  • 购买商品未入库
  • 电脑开机后无显示,但主机电源指示灯长亮
  • 主营业务成本和库存商品区别
  • 质押的应收票据怎么做账
  • 税盘维护费的账务处理
  • php数组有哪几种类型
  • php 抓取别的网站的内容
  • kms.exe
  • PHP:pcntl_sigprocmask()的用法_PCNTL函数
  • php数据迁移
  • 医院产生的相关法律法规
  • php echo js
  • 本月职工工资
  • 前端axios是什么
  • HTML+CSS+JS+Jquery+练手项目+...合集(前端学习必备,持续更新中...)
  • golang 和 java
  • 公司租赁个人车辆怎么开发票
  • 应交税费-应交增值税
  • 利息应怎么录入收入
  • 固定资产减值损失计入
  • 税控系统技术维护费全额抵扣分录
  • 相同的商品附带不同的赠品发布
  • 防伪税控技术维护费普通发票怎么申报
  • 捆绑销售如何做会计处理合适?
  • 施工企业预估成本怎么算
  • 更衣柜属于什么费用
  • 2020年扶贫拨款
  • 编制现金流量表应以什么为基础
  • 多缴所得税返还会计分录
  • 盘亏机器设备
  • 月末结转后应交税费应交增值税一般无余额
  • 母公司销售给控股子公司
  • 采购入库单如何弃审U8
  • 资产负债表本期没有发生额怎么填
  • 预付账款怎么转
  • windows10已经阻止此软件
  • win10 rs3
  • windows 进度条
  • qtaet2s.exe - qtaet2s是什么进程 有什么用
  • windows中常用的菜单有哪三个
  • xp系统如何设置
  • npscheck.exe - npscheck是什么进程 有什么用
  • ssh permission denied password
  • win8怎么切换界面
  • 升级win10系统错误代码0x80072F8F
  • Apache 2.0.55 for Linux 下载
  • win安装ie8
  • jquery.js
  • django orm外键
  • jquery密码验证
  • 微信公众号摇号软件
  • javascript中有哪些数据类型
  • html标签页效果
  • unity二段跳
  • js实现的简洁二次函数
  • 辽宁省国家税务局电子税务局官网
  • 电子发票和普通发票哪个好
  • 国家实行什么制度鼓励电力用户合理调整用电负荷
  • 国地税发展历程
  • 上海买新房办贷款流程
  • 纪检组长如何监督党员
  • 免责声明:网站部分图片文字素材来源于网络,如有侵权,请及时告知,我们会第一时间删除,谢谢! 邮箱:opceo@qq.com

    鄂ICP备2023003026号

    网站地图: 企业信息 工商信息 财税知识 网络常识 编程技术

    友情链接: 武汉网站建设