信号处理中的测量精度与频谱细化问题及其仿真-博文对应的代码

%% 对FFT & 补0后FFT 两种处理方法进行对比 %% %1、画出时域数据 %2、给出两种处理方法下的频域结果 %3、对于单频点，加20dB高斯噪声，给出1000次蒙特卡罗方法下两种处理方式的测量准确度和测量精度的对比% %author: https://blog.csdn.net/xhblair?spm=1010.2135.3001.5343% clear; close all; clc; %% %--------------------------看看FFT&补0 FFT两种处理结果下的频域对比---------------------------% %% % %------------产生时域数据---------------% ftest = [40 44 100 101.5 160 161]; %想要分析的信号中的频点 fs = 400; %采样率设置 N = 256; %采样点数 T = N/fs; %采样的总时间 df = 1/T; %频率分辨率 fprintf(['df = ',num2str(df),'\n']); t = linspace(0,T,N); %设置采样时间 这里的设置应该与后面的频率坐标设置相吻合 % t = (0:1/fs:(N-1)/fs); %设置采样时间（这种设置和上面的设置是不一样的！） x = zeros(1,N); for ii = 1:length(ftest) x = x + sin(2*pi*ftest(ii)*t); end x = awgn(x,20); %加噪声 这里是加高斯白噪 %看看回波数据 figure(1); plot(t,x);xlabel('time/s');ylabel('幅值');title('时域数据');grid on; %---------直接FFT&补0FFT处理----------% fftOut_tmp1 = fft(x,N); fftOut1 = fftOut_tmp1(1:N/2)./N*2; %补0倍数(直接在FFT中增加fft的点数，对应自动在时域的末尾补0) D = 5; fftOut_tmp2 = fft(x,N*D); fftOut2 = fftOut_tmp2(1:N*D/2)./N*2; %看看FFT的结果 xaxis1 = linspace(0,fs/2,N/2); xaxis2 = linspace(0,fs/2,N*D/2); figure(2); plot(xaxis1,abs(fftOut1));xlabel('频率');ylabel('幅值');title('FFT&补0后FFT频域结果对比'); hold on; plot(xaxis2,abs(fftOut2));hold off;grid on;legend('直接FFT','补充4倍数据长度的0后的FFT'); %} %% %----------------针对单频加噪信号进行蒙特卡罗实验，并获取两种方法下的测量准确度和精度---------------------% %% %-----重新构造单频点下的信号-----% ftest = 100; %单频点信号，也是真实频点。 fs = 400; %采样率设置 N = 256; %采样点数 T = N/fs; %采样的总时间 t = linspace(0,T,N); %这里用linspace设置，后面应该也要用linspace设置 x = sin(2*pi*ftest*t); %看看回波数据 figure(5); plot(t,x);xlabel('time/s');ylabel('幅值');title('时域数据');grid on; %-----再在外围增加一个不同SNR下的测量精度的对比-----% SNR = [5 10 15 20 25 30 35 40]; %--------蒙特卡罗仿真----------% 用均方根误差来衡量测量精度 nums_mtkl = 1; %次数 fmeasure_FFT = zeros(length(SNR),nums_mtkl); %装载直接FFT下每次的测量值 fmeasure_Add0FFT = zeros(length(SNR),nums_mtkl); %装载补充0FFT下每次的测量值 D = 5; %补零的倍数 a1 = linspace(0,fs/2,N/2); df_fft = a1(2)-a1(1); a2 = linspace(0,fs/2,N*D/2); df_Add0fft = a2(2)-a2(1); accuracy_FFT = zeros(1,length(SNR)); %装载直接FFT下每个SNR值下的误差值(测量准确度)。 accuracy_Add0FFT = zeros(1,length(SNR)); precision_FFT = zeros(1,length(SNR)); %装载直接FFT下每个SNR值下的均方误差值(测量精确度)。 precision_Add0FFT = zeros(1,length(SNR)); for kk = 1:length(SNR) tmp1 = 0; %用来装在直接FFT下计算均方根误差的中间值 tmp2 = 0; for ii = 1:nums_mtkl x_use = awgn(x,SNR(kk)); %加噪声 % x_use = x + (1 + SNR(kk)*randn(1,length(x))); fftOut_tmp1 = fft(x_use,N); fftOut1 = abs(fftOut_tmp1(1:N/2))./N*2; index1 = find(fftOut1 == max(fftOut1)); %因为只有一个频点，所以可以通过直接找到最大值的方式来得到测量频点。但是每次似乎都是同一个点？？ % fmeasure_FFT(kk,ii) = index1*fs/N; %这里频率的计算是不是有问题？步进如果要对应前面信号的设置，步进应该不是这样算的。 fmeasure_FFT(kk,ii) = (index1-1)*df_fft; tmp1 = tmp1 + (fmeasure_FFT(kk,ii) - ftest)^2; fftOut_tmp2 = fft(x_use,N*D); fftOut2 = abs(fftOut_tmp2(1:N*D/2))./N*2; index2 = find(fftOut2 == max(fftOut2)); %因为只有一个频点，所以可以通过直接找到最大值的方式来得到测量频点。 % fmeasure_Add0FFT(kk,ii) = index2*fs/N/D; fmeasure_Add0FFT(kk,ii) = (index2-1)*df_Add0fft; tmp2 = tmp2 + (fmeasure_Add0FFT(kk,ii) - ftest)^2; % figure(30); % % subplot(122); % xaxis1 = linspace(0,fs/2,N/2); xaxis2 = linspace(0,fs/2,N*D/2); % plot(xaxis1,fftOut1); hold on; plot(xaxis2,fftOut2);hold off; xlabel('频率');ylabel('幅值'); % legend('直接FFT','补充4倍数据长度的0后的FFT');title('FFT&补0后FFT频域结果对比'); grid on; % subplot(121) % plot(fftOut1); hold on; plot(fftOut2);hold off; xlabel('索引');ylabel('幅值'); end %分别计算误差以及均方根误差% accuracy_FFT(kk) = mean(fmeasure_FFT(kk,:)) - ftest; accuracy_Add0FFT(kk) = mean(fmeasure_Add0FFT(kk,:)) - ftest; precision_FFT(kk) = sqrt(tmp1/nums_mtkl); precision_Add0FFT(kk) = sqrt(tmp2/nums_mtkl); end %计算标准差 % datatmp1 = 0; % datatmp2 = 0; % for jj = 1:nums_mtkl % datatmp1 = datatmp1 + (fmeasure_FFT(jj) - mean(fmeasure_FFT))^2; % datatmp2 = datatmp2 + (fmeasure_Add0FFT(jj) - mean(fmeasure_Add0FFT))^2; % end % precision_FFT = sqrt(datatmp1/nums_mtkl); % precision_Add0FFT = sqrt(datatmp2/nums_mtkl); %因为每次测量的结果是一样的，所以每个SNR下的均值误差其实等于标准差... %画图% figure(3); subplot(221); plot(SNR,accuracy_FFT);xlabel('SNR');ylabel('均值误差(测量准确度)');title('不同SNR下，直接FFT处理时的测量准确度'); subplot(222); plot(SNR,precision_FFT);xlabel('SNR');ylabel('均方根值(测量精度)');title('不同SNR下，直接FFT处理时的测量精度度'); subplot(223); plot(SNR,accuracy_Add0FFT);xlabel('SNR');ylabel('均值误差(测量准确度)');title('不同SNR下，补充0后FFT处理时的测量准确度'); subplot(224); plot(SNR,precision_Add0FFT);xlabel('SNR');ylabel('均方根值(测量精度)');title('不同SNR下，补充0后FFT处理时的测量精度');