数字信号处理频域采样及恢复（离散频谱到连续频谱） MATLAB

Posted 2023-05-01

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了数字信号处理频域采样及恢复（离散频谱到连续频谱） MATLAB相关的知识，希望对你有一定的参考价值。

请教用MATLAB 如何把满足频域采样定理的X（K）恢复到X（ejW）　　知道频域内插公式　　但是不会写　请大神帮忙　最好给出源程序

参考技术A clc
clear;
close;
fs = 8000;%采样频率
N = 256;%采样点数
T = fs/N;%频率分辨率
%deltf = fs/N = 8000/256 = 31.25
%输入信号
f1 = 1000;
t = 0:1/fs:(N-1)/fs;
x = 0.2*sin(2*pi*t*f1);
w = 0:fs/N:(fs-(fs/N));
%计算DFT
rol = exp(-1i*2*pi/N); %旋转因子
dft = zeros(1,N);
sum = 0;
for k = 1:N
for i =0:N-1
dft(k) = dft(k) + x(i+1)*rol^(k*i);
end
end
%频域插值提高分辨率
w1 = 0:0.5*0.001*2*pi:2*pi-0.5*0.001*2*pi;
dtft = zeros(1,length(w1));
asum = 0;
%插值程序
for i =0:length(dtft)-1
w = 2*pi*i/length(dtft);
for k = 0:N-1
asum = dft(k+1)*fik(w,k,N);
dtft(i+1) = dtft(i+1)+asum;
end
end
w2 = 0+T:0.5*fs*0.001:fs-0.5*fs*0.001+T;
w3 = 0:fs/N:(fs-(fs/N));

%2000点fft

real = fft(x,2000);
w4 = 0:fs/2000:fs-4;

%频域插值提高分辨率
% w1 = 0:0.5*0.001*2*pi:2*pi-0.5*0.001*2*pi;
% dtft = zeros(1,length(w1));
% asum = 0;
% %插值程序
% for i =0:length(dtft)-1
% w = 2*pi*i/length(dtft);
% for k = 0:N-1
% asum = dft(k+1)*fik(w,k,N);
% dtft(i+1) = dtft(i+1)+asum;
% end
% end
w2 = 0+T:0.5*fs*0.001:fs-0.5*fs*0.001+T;
w3 = 0:fs/N:(fs-(fs/N));
%加窗插值处理
nfft = 2000;
win = hanning(length(x));
x_hannig = x.*win';
y2 = fft(x_hannig,N);
y3 = fft(x_hannig,2000);
w1 = 0:0.5*0.001*2*pi:2*pi-0.5*0.001*2*pi;
dtft = zeros(1,length(w1));
asum = 0;
%插值程序
for i =0:length(dtft)-1
w = 2*pi*i/length(dtft);
for k = 0:N-1
asum = y2(k+1)*fik(w,k,N);
dtft(i+1) = dtft(i+1)+asum;
end
end

w5 = 0:0.5*fs*0.001:fs-0.5*fs*0.001;
figure
plot(w5,abs(dtft),'o',w4,abs(y3),'g')
legend('插值图像','理想逼近图像')
xlabel('频率')
ylabel('幅值')
title('fs = 8000hz fc=1000hz加汉宁窗结果') 参考技术B 我也在问嘤嘤嘤

频域分析频谱泄露频率分辨率栅栏效应

一、时域加窗

现实生活中的信号大部分是连续的，通过对连续的信号进行采样得到散时间信号，但是计算机所能处理的数据都是有限长的，因而我们可以对原始序列做加窗处理使其成为有限长序列。
以矩形窗为例，其时域表达式为：
式(1)中， $N = M + 1$ ，为矩形窗的长度。

对无限长序列进行加窗处理，就是对序列在时域上乘以一个窗函数。
由卷积定理可以得到，时域的相乘等于频域的卷积。

设仿真信号的时域表达式为：
$x(t)=A_{0}*cos(2πf_{0}t)+A_{1}*cos(2πf_{1}t)$
$x (t)$ 做傅里叶变换（FT）的频域表达式为：
$X(jΩ)=A_{0}πδ(Ω+Ω_{0})+A_{0}πδ(Ω-Ω_{0})+A_{1}πδ(Ω+Ω_{0})+A_{1}πδ(Ω-Ω_{0})$
连续信号 $x (t)$ 的波形及频谱如图1所示。

连续信号 $x (t)$ 经过采样后，得到的离散时间的表达式为：
$x[n]=x(t)|_{t=nT_{s}}$
离散序列 $x [n]$ 做离散时间傅里叶变换（DTFT）的频域表达式为：
$X(e^{jw})=\\frac{1}{T_{s}}\\sum_{k=-∞}^{∞}X(j\\frac{w}{T_{s}}-jk\\frac{2π}{T_{s}})$
离散序列 $x [n]$ 的波形及频谱如图2所示。

矩形窗函数 $w [n]$ 做离散时间傅里叶变换（DTFT）的频域表达式为：
$W(e^{jw})=e^{-jw(N-1)/2} *\\frac{sin(wN/2)}{sin(w/2)}$

矩形窗函数 $w [n]$ 的波形及频谱如图3所示。

离散序列 $x [n]$ 与窗函数 $w [n]$ 的卷积为：
$V(e^{jw})=\\frac{1}{2π}\\int_{-π}^{π}X(e^{jθ})W(e^{j(w-θ)})dθ=\\frac{A_{0}}{2}W(e^{j(w+w_{0})})+\\frac{A_{0}}{2}W(e^{j(w-w_{0})})+\\frac{A_{1}}{2}W(e^{j(w+w_{0})})+\\frac{A_{1}}{2}W(e^{j(w-w_{0})})$
截断后的离散序列 $v [n]$ 的波形及频谱如图4所示。