[M] 修改栅栏效应的测试方案：频率分辨率测试

1. 混叠测试错别字：Dedimate --> Decimate； 2. 系统辨识测试将多次仿真汇总到同一个图； 3. 栅栏效应将单频做谱改为两个频率做谱，有利于对比频率分辨率。
iChunyu · Oct 22, 2020 · 82bdffa · 82bdffa
1 parent 19b332f
commit 82bdffa
Show file tree

Hide file tree

Showing 3 changed files with 29 additions and 27 deletions.
diff --git a/Demo03_SysIdentify.m b/Demo03_SysIdentify.m
@@ -17,23 +17,20 @@
 fs = 10;
 T = 5e3;
 N = T*fs;
+t = (0:N-1)'/fs;
 
+figure('Name','Identify')
 for k = 1:5
-    t = (0:N-1)'/fs;
-    xn = randn(N,1)*sqrt(fs/2);
-
+    xn = randn(N,1)*sqrt(fs/2); 
     y = lsim(sys,xn,t);
     [pxx,f] = iLPSD(y,fs);
-
-    [A,~,~] = bode(sys,2*pi*f);
-    A = squeeze(A);
-
-    figure('Name','Identify')
-    loglog(f,sqrt(pxx),'DisplayName','Identified')
+    loglog(f,sqrt(pxx),'DisplayName',sprintf('Identified%0d',k))
     grid on
     hold on
-    loglog(f,A,'DisplayName','Real')
-    legend
-    xlabel('Frequency (Hz)')
 end
+[A,~,~] = bode(sys,2*pi*f);
+A = squeeze(A);
+loglog(f,A,'Color','k','LineStyle','--','DisplayName','Real')
+legend
+xlabel('Frequency (Hz)')
 ylabel('Gain')
diff --git a/Test02_Aliasing.m b/Test02_Aliasing.m
@@ -34,19 +34,19 @@
 %% 10倍直接抽取，采样率降为100Hz
 x1 = x0(1:10:end);
 [pxx,f] = iLPSD(x1,fs/10);
-loglog(f,sqrt(pxx),'DisplayName','Dedimate-100Hz')
+loglog(f,sqrt(pxx),'DisplayName','Decimate-100Hz')
 legend('location','southwest')
 
 %% 50倍直接抽取，采样率降为20Hz
 x2 = x0(1:50:end);
 [pxx,f] = iLPSD(x2,fs/50);
-loglog(f,sqrt(pxx),'DisplayName','Dedimate-20Hz')
+loglog(f,sqrt(pxx),'DisplayName','Decimate-20Hz')
 legend('location','southwest')
 
 %% 100倍直接抽取，采样率为10Hz
 x3 = x0(1:100:end);
 [pxx,f] = iLPSD(x3,fs/100);
-loglog(f,sqrt(pxx),'DisplayName','Dedimate-10Hz')
+loglog(f,sqrt(pxx),'DisplayName','Decimate-10Hz')
 legend('location','southwest')
 
 %% 抗混叠滤波后抽取，采样率为10Hz

diff --git a/Test03_FenceEffect.m b/Test03_FenceEffect.m
@@ -14,49 +14,54 @@
 
 %% 构造单频信号
 fs = 100;
-fsig = 1;
+fsig1 = 1;
+fsig2 = 1.001;
 T = 1000;
 N = T*fs;
 
 t = (0:N-1)'/fs;
-x = sin(2*pi*fsig*t);
+x1 = sin(2*pi*fsig1*t);
+x2 = sin(2*pi*fsig2*t);
+x = x1+x2;
 
 figure('Name','Sine')
-plot(t,x)
+plot(t,x1)
+hold on
 grid on
+plot(t,x2)
 xlabel('Time (s)')
 ylabel('Value')
-xlim([0 5])
+xlim([0 10])
 
 %% 不同频率点计算PSD
 win = rectwin(N);
-df = 0.005;
+df = 0.003;
 
 figure('Name','PSD')
 [pxx,f] = periodogram(x,win,N*100,fs,'onesided');
 plot(f,sqrt(pxx),'DisplayName','DTFT')
 hold on
 grid on
 legend
-xlim([fsig-df fsig+df])
+xlim([fsig1-df fsig2+df])
 xlabel('Frequency (Hz)')
 ylabel('Power Spectrum')
 
 %%
-markersize = 30;
-[pxx,f] = periodogram(x,win,N*10,fs,'onesided');
-plot(f,sqrt(pxx),'DisplayName','10nfft',...
+markersize = 25;
+[pxx,f] = periodogram(x,win,N/2,fs,'onesided');
+plot(f,sqrt(pxx),'DisplayName','nfft/2',...
     'Marker','.','MarkerSize',markersize,'LineStyle','none')
 legend
 
 %%
-[pxx,f] = periodogram(x,win,5*N,fs,'onesided');
+[pxx,f] = periodogram(x,win,N,fs,'onesided');
 plot(f,sqrt(pxx),'DisplayName','nfft',...
     'Marker','.','MarkerSize',markersize,'LineStyle','none')
 legend
 
 %%
-[pxx,f] = periodogram(x,win,N/2,fs,'onesided');
-plot(f,sqrt(pxx),'DisplayName','nfft/2',...
+[pxx,f] = periodogram(x,win,N*5,fs,'onesided');
+plot(f,sqrt(pxx),'DisplayName','5nfft',...
     'Marker','.','MarkerSize',markersize,'LineStyle','none')
 legend