互功率谱
类似从自相关函数推出互相关函数,也可以从单个随机过程的功率谱推出互功率谱。两个随机过程X(t)和Y(t),它们的互功率谱定义是
如果X(t)和Y(t)是联合平稳的,则
绝对可积。则它们的互功率谱和互相关函数是一对傅里叶变换对:
互功率谱的性质:
非平稳过程的功率谱
广义功率谱密度
设非平稳随机过程X(t)的自相关函数RX(t1,t2)是两个时间点的函数,则
的二维傅里叶变换
为X(t)的广义功率谱密度,逆变换为
广义谱的定义在数学上与平稳随机过程是相通的。假设
仅在两者相等处有值,两者不等时的结果是0,则
这就是平稳情况下的功率谱。
时变功率谱
假设t1=t,t2=t+τ,则非平稳过程的自相关函数
是t和
的函数,时变谱定义成
如果采用对称相关函数
则这种时变谱称为维格纳-威利谱(WV谱)其中,
还可以写成
其中
称为维格纳分布,非平稳随机过程X(t)的WV时变谱是该过程维格纳分布的数学期望。
例子
设噪声调制的振荡信号
其中N(t)是平稳噪声,求X(t)的功率谱。
观察这个信号,发现不同时刻的信号(可能的)均值受到了三角函数的调制。X(t)是非平稳随机过程。求功率谱的方法一般就是求相关函数的傅里叶变换,而对于非平稳过程,有广义功率谱密度、时变功率谱等表示形式,但最常用的还是对即时相关函数做时间平均,然后进行傅里叶变换,如下。
所谓的对时间函数求平均,指的是
对上述自相关函数求时间平均就得到了结果:
matlab互相关函数估计
参见文档https://www.mathworks.com/help/matlab/ref/xcorr.html
>> a=0.8;
>> sigma=2;
>> N=500;
>> u=randn(N,1);
>> x(1)=sigma*u(1)/sqrt(1-a^2);
>> for i=2:N
x(i)=a*x(i-1)+sigma*u(i);
end
>> R=xcorr(x,'coeff');
>> plot(R,'r','linewidth',2);
matlab功率谱估计
参照函数https://www.mathworks.com/help/signal/ref/periodogram.html
例如,估计两个正弦信号加正态白噪声的功率谱,信号
其中,
>> t=0:0.001:0.3;
x=cos(2*pi*t*300)+cos(2*pi*t*310)+randn(size(t));
subplot(2,2,1);
periodogram(x,[],512,10000);
axis([0 500 -50 0]);
subplot(2,2,2);
window = hann(301);
periodogram(x,window,512,1000);
axis([0 500 -50 0]);
>> R=xcorr(x)/150000;
>> Pw=fft(R);
>> subplot(2,2,3);
>> f=(0:length(Pw)-1)*1000/length(Pw);
>> plot(f,10*log10(abs(Pw));
>>> plot(f,10*log10(abs(Pw)));;
>> axis([0 100 -50 0]);
>> grid on
>> subplot(2,2,4);
>> pwelvh(x,128,64,[],1000);
>>> pwelch(x,128,64,[],1000);
>> axis([0 500 -50 0]);
不加窗的周期图法如下图。就是先对离散序列做傅里叶变换,再将功率谱估计成频域能量的时间平均。
则功率谱估计成是
实际情况中数据截取的长度是有限的,数据的截断会导致一定截断误差的产生。减少截断效应,常采用数据加窗。下图是加汉宁窗后的周期图功率谱估计。
也可以采用传统的自相关函数法估计(先求出相关函数的估计,再对相关函数做傅里叶变换)
最后给出matlab另一个自带的韦尔奇功率谱估计。
例子:数字图像直方图
随机过程可以随时间变化,也可以看成随空间变化。数字图像可以看作随位置变化的随机序列。
数字图像灰度级的直方图,指的是反映图像中灰度级与出现这种灰度之间概率的图形。设R代表图像中像素灰度级,R的取值范围是[0,L-1],L为总的灰度级数,
:第k个灰度级,越往下越黑。nk是具有灰度级
的像素度。直接想到可以对灰度做归一化,
f(rk)是灰度级出现的概率估计。归一后就有了
直方图可用于图像压缩、图像增强等处理技术中。下面是最最简单的处理:
>> I=imread('微信图片_20200427230532.jpg');
>> figure
subplot(1,2,1)
imshow(I)
subplot(1,2,2)
imhist(I,64)
调亮:I=I+30;
直方图灰度级集中向高端移动。
均衡(从上面图的像素灰度级分布转换成服从均匀分布的灰度级。)
不妨先假设R是连续的,它(灰度级)归一化了,概率密度
,对R做变换得到新变量S
转移函数是0到1之间的单调函数。S的概率密度
比较经典的是积分一下,相当于求分布函数
可以证明S在0到1服从均匀分布(利用随机变量的分布函数求解其反函数可得到任意分布随机数,随机变量的函数变换可获得任意分布随机数),
。
这里把R变成离散的之后,就是把上面概率密度变成概率,把积分改成求和。R=rk的概率是
那同理,
这样,图像像素由R变成S,就服从均匀分布了。matlab仍然已经写好了这个过程。
>> J=histeq(I);
>> figure
subplot(1,2,1)
imshow(J)
subplot(1,2,2)
imhist(J,64)
随机过程应用:判断说话中的清音浊音
衡量标准综合如下:
语音片段的能量:短时能量通过大约长度为20ms的片段进行衡量。显然,浊音片段往往比清音片段有更高的能量,这从波形中就很容易看出。
过零率(zero crossing rate. ZCR):
在昨天的例子中实际上已经实现了:
#在导入一段音频后,取一小段分析
n0= 57000
n1 = 57100
plt.figure(figsize = (14,5))
plt.plot(x[n0:n1])
plt.grid()
计算过零点个数:
zero_crossings = librosa.zero_crossings(x[n0:n1],pad = False)
print(sum(zero_crossings))
往往,浊音部分的过零率小 ,清音部分的过零率大(频率往往比较高)。而沉静片段的过零率往往和说话的过零率在一个量级,都比清音的低。当然,沉静片段的能量会小很多,所以可以区别开来。这里举一个例子,我说“东风来了,春天的脚步近了”,然后让它判断清音和浊音大概怎么分类。
首先读入文件并画图。读取前8000个样本点。
>> a=fread(fp,8000);
>> plot(a)
>> fp=fopen('春天的脚步近了.m4a');
fseek(fp,60244,-1);
>> a=fread(fp,8000);
>> plot(a)
像刚才在python中做的一样,计算零交叉点的数目。8000确实有点多了,这里不妨取第5000到第7000的样本点进行处理。(已事先在音频软件中观测到5000到7000的样本是在清音阶段)
figure;
for i=5000:5199
if (a(i)>128) && (a(i+1)<128)
x=x+1;
else x=x;
end
end
disp('number of zero crossing='); disp(x);
再在音频软件中定位到春天an处,(浊音)作出这附近的曲线:
fseek(fp,65042,-1);
b=fread(fp,2000);
plot(b);
用相同的方法,得到零交叉点数目,可见大概只有上个数目的不到一半。
再在音频软件中定位到脚步的j处(清音),作出这附近的曲线:
fseek(fp,6724,-1);
c=fread(fp,200);
plot(c);
可以看出,它的振幅变化幅度较浊音更大。
归一化自协方差
这个量规定了相邻音乐/语音片段的相关性。讲人话,C1就可以理解成这时求的相邻单位时间延迟的自协方差。对于浊音,有大量低频信号能量,在浊音附近是的信号高度自相关的,这个值很大(讲人话,低频信号说明信号的乱起八糟的高频谐波占比比较少,信号比较“单纯”,相互的依赖性比较强);而在清音区域自协方差就很小,在寂静区域会更小(讲人话,就是噪声基本上无法预测、是独立产生的)。与之类似,在语音信号中通常再引入一个叫谱倾斜的东西,定义成
可以用matlab对“春天来了”做差分,求出它的谱倾斜度(spectrum tilt)(当然,需要事先录制"春天来了"4字)
>> fp=fopen('录音 (5).m4a');
%fseek(fp,800,-1);
a=fread(fp,400,'short');
a=a-128;
subplot(2,1,1);plot(a);
xlabel('sample no.');ylabel('amplitude');
for j=1:200,
fseek(fp,4*j,-1);
a=fread(fp,100);
a=a-128;
sum(j)=0;
for i=2:100,
sum(j)=sum(j)+(a(i)*a(i-1));
end
sum1(j)=0;
for i=2:100,
if a(i)==0,
a(i)=0.1;
end
sum1(j)=sum1(j)+a(i)*a(i);
end
s(j)=sum(j)/sum1(j);
end
subplot(2,1,2);plot(s);
得到的是下图。
可以看出,由于我说的比较波澜不惊,波形很规整的了。音频信号的自相关系数(谱倾斜度和自相关函数实际上差不多是同一种量度)在浊音处是十分高的;由于春和天是连在一起的,所以得到的相关性都比较大(接近1);而后来说来字的时候,音调和共振峰都发生了比较强烈的改变,所以它的相关性急剧下降。