阿英讲频率f，角频率Ω和数字频率w的物理含义附MATLAB仿真

转自：http://anony3721.blog.163.com/blog/static/51197420111129503233/

 

古人云：基础不牢，地动山摇。勿在浮沙筑高台。此话真不假，比如MATLAB中下标从1开始而物理概念t从0开始，结果往往会差一点，做FFT后结果会莫名其妙的差一点，做仿真的时候经常会因为这样一些基本概念不清而导致对结果无法正确的解释。盲目的追求多学习，不求甚解是得不偿失的，最后无知的还是你自己。一定要动脑子想想，把知识消化了才能灵活运用。本文是数字信号处理的基本功，是本人学习思考后的总结，网上没有发现有人讲此很基本的内容，相信肯定有不懂的，所以贴出来希望大家能受益。最后，原大家得大智慧，断贪嗔痴，阿弥陀佛。

1。 模拟角频率Ω：单位rad/s大OMEGA的物理含义是2*pi的时间段里面包含y=sin(OMEGA*t)正弦信号的周期数。Ω = 2*pi/T ，现象下手指头绕原点做圆周运动，经过一个周期T的时间，我们绕了Ω 圈，则在2*pi的时间段内，正弦y=sin(OMEGA*t)就会有Ω个完整的波形。我们往往看到Ω = 2*pi/T每秒经历多少弧度，单位rad/s，你想想到Ω 当初其实就是刻画你绕原点画圆圈的快慢了吗？正弦信号和余弦信号实际上是绕圆周运动的点在x轴和y轴上的投影。Ω=2*pi*f

 

t=0:pi/50:2*pi;
for OMEGA = 1:4
y(:,OMEGA) = sin(OMEGA*t).';
str{OMEGA}=['OMEGA=',num2str(OMEGA)];
end
h=plot(t',y);
legend(h,str);

 

 2. 频率f：单位Hz，频率f的物理含义是1s的时间段内包含有f个y=sin(f*t)完整周期的信号波形。根据f = 1/T 可以看出 f 表征的是1s时间段内振动了多少次。频率f不同于角频率Ω是绕着一个周期T时间间隔内物体绕原点转的圈数。由Ω = 2*pi*f ，看出f已经和圆周2*pi脱离关系了，手指头绕原点画圆圈投影到x轴或者y轴上的情况是指尖左右移动，或者上下振动的情况，这就是f 的物理含义：在周期T内上下振动的次数，或者左右移动的次数。T = 1/f 下图中是正弦的情况，f = 1 就是1s内指尖从(1, 0)点向(-1,0)点移动再移回来1次。。。 f = 4 就是1s内指尖从(1, 0)点向(-1,0)点移动再移回来4次

clear;
t=0:1/100:1;
for f= 1:4
y(:,f) = sin(2*pi*f*t).';
str{f}=['f=',num2str(f)];
end
h=plot(t',y);
legend(h,str);

 3. 数字频率w：单位是rad。学习数字信号处理的要明白，数字频率实际上是和采样周期Ts联系在一起的，离开采样周期Ts或者采样频率Fs = 1/Ts单独谈数字频率w是没有实际价值的，因为此时它没有任何物理意义，这个我当初学的时候也没有弄明白，此点大家要谨记！数字频率w是从单位圆上的N点等间隔采样而来的，这个N不是别的就是数字周期，先给出数字频率w和数字周期N之间的关系：N=(2*pi/w)*k当信号是 以N为周期的时候，要求(2*pi/w)是个有理数，对这个不明白的请看数字信号处理的课本。w=2*pi*k/N，可见w的物理含义是相邻的两个采样点之间的弧度，w = ΩT = Ω/Fs 是用Fs归一化后的频率。此时想想手指不沿着单位圆绕着原点做转圈了，而是在单位圆上有N个树坑，第一个树坑和第二个树坑之间的弧度大小为w。注意：采样频率Ωs >= 2*Ω 和数字频率是一点关系都没有！（此处Ωs=2*pi*Fs ）顺便说下数字周期N，N有两种含义：①f(n) = f(n+N)

 ② 周期信号 f(n)的一个周期内有0~N-1共N个采样点。 cos(2pi*f*t) = cos(Ω *t) = cos(Ω*n*Ts) = cos(Ω*Ts*n) = cos(w*n)

clear;
w = pi/4;
N = 2*pi/w; % N = 8 ，此处是在数据上采样N个点，FFT时在频谱上也采样N个点
n = 0:N-1;
x = sin(n*w);
h=plot(n, x, '-o');

 注意上图中没有时间信息，所以没有任何物理含义。这在实际中是没有用的，所以要加上采样信息Ts或者Fs，才能知道具体的物理含义,比如这些圆周上的树坑是以多快的速度挖好的？Ω = w*Fs=w/Ts是1天挖好(Ts较小)的还是3天挖好的（Ts较大）。由于采样丢失了Ts信息，即丢失了挖树坑的快慢信息，而仅仅将单位圆上树坑的个数N以及树坑之间的间距w保留了下来，仅仅是个最终结果。现实世界中有果比有因，非因果的东西是没有意义的，所以我们必须弄清楚这些树坑是以什么多长时间挖好的，即要补充Ts采样信息，才可以对应到模拟角频率Ω，Ω*Ts = w 或者频率 f 上，2*pi*f*Ts = w 。

clear; close all;

% 该信号的数字周期N = 8，模拟周期T = N*Ts = 0.008s ，实际频率 f = 125 Hz

w = pi/4;
N = 2*pi/w; % N = 8
n = 0:N-1; % n = 0:N 可能会更好看一些，但是要清楚第N+1点可是下一个采样周期的第一个点
x = sin(n*w);
h=plot(n, x, '-o'); %注意n没有定标，没有物理含义！

Fs = 1000; %采样频率为1000Hz

Ts = 1/Fs;

t = n*Ts; % t 时间序列的给法永远只有这么一种，请铭记！

T = N*Ts; %模拟周期T

f = w*Fs/(2*pi); % 信号的真实频率f

figure; plot(t, x)

% （1） 做FFT的点数和时间采样的点数N（注意和数字周期N相同）相同时频率f的定标 --- 数据上采样N个点，频谱上采样N个点

freq = n*Fs/N - Fs/2; %频率序列的定标，注意此式子是本来面目，不要给此式子穿个马甲就不认识了，下面是此式子的变形

% freq = (n/N - 1/2)*Fs;        freq = (n - N/2)/N * Fs; % 注意0~N-1时 第N 点对应着 Fs， N/2 点对应着 Fs/2，Fs/N 就是频率分辨率

% freq = (n - N/2)/N;  就是归一化频率，从-0.5 到 0.5 - 1/N , 看看不是 -0.5 ~ 0.5 哦，实际结合硬件编程就要这样干

X = fftshift(abs(fft(x))); %采用 N= 8点的FFT时的FFT，做FFT的点数和时间采样的点数相同

x_IFFT = ifft(fft(x));

figure; plot(freq, X);grid

figure; plot(t, x_IFFT); grid

t_ = (0:N-1)*Ts - N*Ts/2; % 等价的有 t_ = (0-N/2:N-1-N/2)*Ts;  即 t_ = (-N/2: N/2 - 1)*Ts

figure; plot(t_, x_IFFT);grid;

 % 可以看出来频谱没有混叠！而且频率正确

% （2） 做FFT的点数和时间采样的点数N不同时频率f的定标，补零后FFT

补零后的FFT和补零前的FFT，两者会有较大的不同，但两频谱的包络还是一致的。设补零前数据长N，补零后数据长M（补了M-N个零值），则补零前的FFT有N条谱线，分别代表的频率点是（0,1,...,N-1)*fs/N；补零后的FFT有M条谱线，分别代表的频率点是（0,1,...,M-1)*fs/M。由于补零前后数据长度不一样，它们的分辨率（分别为df1=fs/N，df2=fs/M）不一样，在频域中谱线所代表的频率也不一样，所以这两个频谱所描述的对象也不相同。这里举一个例子，fs=1000HZ，补零前后数据长度N=500和M=800，对应的df1=2，df2=1.25。补零前的频谱是对应于0,2,4,...,500HZ的频谱，而补零后的频谱是对应于0,1.25,2.5,...,500HZ的频谱，所以两频谱中对应频率不两同，描述当然不同。注意补零前后fs始终没有变！
但是当特殊情况：M=（2^n)*N时，补零后的频谱相当于在补零前的频谱中插入（2^n)-1条谱线。与补零前的频谱中相重合的谱线，它们的幅值和相位完全一致。

Nfft = 1024; % 此处 Nffft/N = 128，所以要对结果的频谱进行128点的抽取，才可以得到正确的频率定标！

X_Nfft = fftshift(abs(fft(x,Nfft)));
freqNormalized = ((0:Nfft-1) - Nfft/2)/Nfft;

freq = freqNormalized(1:128:end)*Fs; % 结果和上面的结果完全相同 -500  -375  -250  -125     0   125   250   375

X_ = X_Nfft(1:128:end);

figure; plot(freq , X_);grid

4. 说说0~N-1和1~N和归一化频率

归一化频率：  -1 ---------- -1/2 ---------- 0 ---------- 1/2 ----------- 1

频率f：           -Fs --------- -Fs/2 -------- 0 ---------- Fs/2 ---------Fs

模拟角频率Ω：-Ωs --------- -Ωs/2 ------ 0 ----------- Ωs/2 -------- Ωs

数值频率w：    -2*pi --------- -pi --------- 0 ---------- pi ------------- 2*pi

0~N-1和1~N都是描述单位圆上N个树坑的分布情况的，是一个正弦周期内应该采样的点的范围。   

① 0~N-1描述的概念是这样的：第一个采样点N=0对应w的w=0，最后一个采样点N-1点对应w = 2*pi*(N-1)/N。此时w=pi 处对应f = Fs/2或Ω = Ωs/2的模拟最高频分量对应 N/2而不是[0 + (N-1)]/2 = (N-1)/2，因为最后一点即N-1点不对应Fs，从图像上可以很清楚的看到这一点，不注意的话在此就要犯错。A/D采样就是按这样0~N-1的方式操作。比如fftshift后往往要调整下频率范围此时要用到中点(N-1)/2

② 1~N描述的也是单位圆上N个树坑的分布情况，只不过第一点N=1的位置在 w = 2*pi/N的位置，最后一点N点对应w = 2*pi*N/N = 2*pi是零频率的位置。由于MATLAB中下标从1开始，所以这样表达有时候比较方便。此时w=pi 处对应f = Fs/2或Ω = Ωs/2的模拟最高频分量对应N/2，而不是[1 + N]/2，因为N=1点采样没有在t=0处。

下面看一个 1：N的例子

clear;

N=100; % 数字周期 N
Fs=2000; % 采样频率，只有确定了采样频率才能从数字信号频率确定实际模拟频率
Ts=1/Fs;
n=1:N; % n = 0:N 可能会更好看一些，但是要清楚第N+1点可是下一个采样周期的第一个点

w=0.02*pi
x=sin(w*n); %数字频率w=0.02*pi
t=n*Ts; % t 时间序列的给法永远只有这么一种，请铭记！
plot(t,x);grid on; %注意此时信号的第t = 0时没有采样点
T=N*Ts

f = w*Fs/(2*pi); % 信号的真实频率f

% 该信号数字周期为2*pi/w=100，模拟周期为T=100*Ts=0.05s，实际频率为w*F/2*pi=20Hz

freq1 = n*Fs/N - Fs/N; % 1~N时 频率序列 f 的取值应该和 0~ N-1时的情况一样  

X = (abs(fft(x)));

figure; plot(freq1, X);grid

freq = (n-1)*Fs/N - Fs/2; %频率序列

X = fftshift(abs(fft(x)));

figure; plot(freq, X);grid