MATLAB功率谱估计

close all;
clc;
clear all;
fs = 1024;
len = 1024;
t = 0:1/fs:(len-1)/fs;
x = 5.3 * sin(2*pi*16*t);
ffty = fft(x, len);
ffty = abs(ffty) / len;
f = 0:fs/len:fs/2;
plot(f, ffty(1:len/2+1))
title('FFT')
figure
py = ffty .^ 2;
py = py * 2;
py = py ./ (fs/len);
plot(f, py(1:len/2+1))
title('PSD')

figure
[pxx, fxx] = periodogram(x, [], [],  fs);
plot(fxx, pxx)
title('periodogram')

figure
cxn=xcorr(x);%计算序列的自相关函数
cxn=cxn(len:end)+[0 cxn(1:len-1)];
CXk=fft(cxn,len);
Pxx=abs(CXk) * 2 / (len * fs);
plot_Pxx=10*log10(Pxx); 
plot(f, plot_Pxx(1:len/2+1))
title('xcorr')

互功率谱估计

clear;
clc;

fs = 1024;
len = 1024;
t = 0:1/fs:(len-1)/fs;
x = 3.7 * sin(2*pi*16*t);
y = 1.9 * sin(2*pi*32*t);

w = blackman(1024, 'periodic');

[p, f] = cpsd(x, y, w, 0, 1024, fs);
plot(f, abs(p))
title('cpsd')
max(abs(p))

z = xcorr(x, y);
cnx = z(len:end) + [0 z(1:len-1)];
w = w';
x1 = cnx .* w;
cxk1 = fft(x1, 1024);
p1=abs(cxk1) * 2 / (1024 * fs);

w1 = w .* w;
u = sum(w1) / 1024;
%p = p1/u;
p = p1;

f = 0:fs/1024:fs/2;
figure
plot(f, p(1:1024/2+1))
title('xcorr')
max(p)

复倒谱

Fs = 100;
t = 0:1/Fs:1.27;
s1 = sin(2*pi*45*t);
s2 = s1 + 0.5*[zeros(1,20) s1(1:108)];
c = cceps(s2);
plot(t,c)
xlabel('Time (s)')
title('Complex cepstrum')

h = fft(s2);
x = angle(h);
n = length(x);
y = unwrap(x);
nh = fix((n+1)/2);

idx = nh+1; 
if length(y) == 1
  idx = 1;
end
nd = round(y(idx)/pi);
y(:) = y(:)' - pi*nd*(0:(n-1))/nh;


logh = log(abs(h)) + sqrt(-1)*y;
y = real(ifft(logh));
figure
plot(t, y);
title('fft ceps')