多速率采样的demo.rar

% 变速率采样的仿真实验 %% 参数设置 A = [1 2 1.5 ]; f = [10 20 60]; phi = [0 pi/3 pi/2]; N = 1000; t = 2; % 观察时间 n = linspace(0,t-1/N,N); % 采样间隔为Ts = t/N,“种树”问题 x = zeros(1,N); for i = 1:3 x =x+ A(i)*cos(2*pi*f(i).*n+phi(i)); end %% 抽取 M = 2; x1 = zeros(1,ceil(N/M),'double');%分配空间 % 抽取 for i = 1:M:N x1((i+M-1)/M) = x(i); end %% 插值 L = 3; x2 = zeros(1,L*N,'double');%分配空间 % 插0值 for i = 1:N x2(1,(i-1)*L+1) = x(i); end %% 分数采样 假设为3/4 或者叫变速率级联 % 中间还需要滤除镜像 M_ = 4; L_ = 3; N_ = L_*N; tmp = zeros(1,N_,'double'); for i = 1:N tmp(1,(i-1)*L_+1) = x(i); end %滤除镜像 fft_tmp = fft(tmp); fft_tmp(N/2+1:N_-N/2) = zeros(1,length(N/2+1:N_-N/2)); tmp = real(ifft(fft_tmp)); %下采样 x3 = zeros(1,ceil(N_/M_),'double'); for i = 1:M_:N_ x3((i+M_-1)/M_) = tmp(i); end %% fft，以及对应的横轴 X = fftshift(fft(x)); X_amp = abs(X)/length(x); w_axis = linspace(-pi,pi-2*pi/length(x),length(x)); %[-pi,pi),必须是偶数个采样点才行 f_axis = w_axis*length(x)/t/(2*pi); % 根据DTFT和DFT的关系推导的公式 X1 = fftshift(fft(x1)); X1_amp = abs(X1)/length(x1); w1_axis = linspace(-pi,pi-2*pi/length(x1),length(x1)); f1_axis = w1_axis*length(x1)/t/(2*pi); X2 = fftshift(fft(x2)); X2_amp = abs(X2)/length(x2); w2_axis = linspace(-pi,pi-2*pi/length(x2),length(x2)); f2_axis = w2_axis*length(x2)/t/(2*pi); X3 = fftshift(fft(x3)); X3_amp = abs(X3)/length(x3); w3_axis = linspace(-pi,pi-2*pi/length(x3),length(x3)); f3_axis = w3_axis*length(x3)/t/(2*pi); %% 绘图 figure(1) subplot(421); plot(n,x);xlabel("t/s") subplot(422); plot(f_axis,X_amp);axis([-250 250 0 1]);xlabel("f/Hz") subplot(423); plot(linspace(0,t-1/length(x1),length(x1)),x1);xlabel("t/s") subplot(424); plot(f1_axis,X1_amp); title("抽取");axis([-250 250 0 1]);xlabel("f/Hz") legend("M = 2") subplot(425); plot(linspace(0,t-1/length(x2),length(x2)),x2);xlabel("t/s") subplot(426); plot(f2_axis,X2_amp);title("插值"); axis([-250*3 250*3 0 1]) % 这个如果将3改成1，将相当于通过了一个理想低通滤波器 xlabel("f/Hz") legend("L = 3") subplot(427); plot(linspace(0,t-1/length(x3),length(x3)),x3);xlabel("t/s") subplot(428); plot(f3_axis,X3_amp);title("分数采样");axis([-250 250 0 1]) xlabel("f/Hz") legend("L_/M_ = 3/4")