信息处理 - Matlab中FFT后信号的能谱 - 吾爱随笔录

Matlab中FFT后信号的能谱

信息处理 matlab fft 傅里叶变换频谱

2022-02-06 12:20:51

我有一个信号向量 x(t) 及其时间向量。我想获得信号的频率表示，特别是 x(t) 的能谱。

有人可以展示一些光来获得能谱吗？能谱和能谱密度，ESD是否产生相同的结果？

任何帮助深表感谢！

 t = 0:.001:.25;
x = 5*sin(2*pi*50*t) + 1.2*sin(2*pi*120*t);
y = x + 2*randn(size(t));
y=y';
Fs = 1000; % Sampling frequency

% Fast Fourier Transform 
    L = length (y); % Window Length of FFT
    nfft = 2^nextpow2(L); % Transform length , Next power of 2 from length of signal
                          %  Exactly same as nfft = pow2(nextpow2(L)) 

% Take fft, padding with zeros so that length(input) is equal to nfft     
    Ydft = fft(y,nfft)/L;

    y_HamWnd = y.*hamming(L);           %hamming window fft
    Ydft_HamWnd = fft(y_HamWnd,nfft)/L;

%FFT is symmetric, throw away second half

    % Take the magnitude of fft of x
%Amplitude Spectrum 01 (Single Sided)
mY = 2*abs(Ydft(1:nfft/2+1)) ;

%Amplitude Spectrum 02 (Single Sided HamWmd)
mY_HamWnd = 2*abs(Ydft_HamWnd(1:nfft/2+1));

%Amplitude Spectrum 03 (Single Sided; Exclude DC-- FoldingFreq)
mY2 = abs(Ydft);
mY2= mY2(1:nfft/2+1);
mY2 (2:end-1) = 2* mY2 (2:end-1); 


    %To get the Power
% P ; single sided; then coefficient square**********Multiply by 2 ?
YdftCoeffi_Sqr = Ydft(1:nfft/2+1).^2;   

% P = y*(conj y) == y^2  ;     %********DO WE NEED TO MUL 2 ?
Pyy_NormNfft = Ydft.*conj(Ydft);   

% P ; single sided; First square, then Multiply by 2
m1YL = abs(Ydft(1:nfft/2+1));  
mSqrYL = m1YL.^2;  
m2YL = mSqrYL*2; 

% Frequency vector
f1 = Fs /2 * linspace(0,1, nfft/2 + 1);
f2 = (0:nfft-1)*(Fs/nfft);    

% Generate the plot, title and labels. 

figure;
plot(t,y); title('Time Domain'); xlabel('Time(s)'); 

figure;
subplot(3,1,1)
plot(f1,mY)
title('Amp Spec 01 --mY = 2*abs(Ydft(1:nfft/2+1)) /L; '); xlabel('Frequency (Hz)'); ylabel('|mY|')
        %axis ([0 500 0 5]); %Zoom in 
subplot(3,1,2)
plot(f1,mY_HamWnd)
title('Amp Spec 02 -- mYHamWnd =2*abs(YdftHamWnd(1:nfft/2+1))/L;'); xlabel('Frequency (Hz)'); ylabel('|mY_HamWnd|')
        %axis ([0 500 0 5]); %Zoom in 
subplot(3,1,3)
plot(f1,mY2)
title('Amp Spec 03 DC-- FoldingFreq -- mY2 ;First, single sided; then square, lastly Multiply by 2'); xlabel('Frequency (Hz)'); ylabel('|mY2|')
        %axis ([0 500 0 5]); %Zoom in 

figure;
subplot(3,1,1)
plot(f1,YdftCoeffi_Sqr); 
title('Ydft Complex Coefficient Square; Single Sided-- YdftCoeffi_Sqr = Ydft(1:nfft/2+1).^2');  xlabel('Frequency (Hz)'); ylabel('YdftCoeffi_Sqr');
subplot(3,1,2)
plot(f1,m2YL); 
title('Ydft Complex Coefficient Magnitude Square; Single sided; then Multiply by 2'); 
xlabel('Frequency (Hz)'); 
ylabel('Power m2YL');
        %axis ([0 500 0 5]); %Zoom in

1个回答

引用的段落给出了一种获得特定频率谐波能量的方法。如果您感兴趣的频率是，则取一段长度为的样本，其中并计算点 FFT。 th FFT bin的幅度来找到二次谐波的能量。要考虑估计的误差，您可以将 bin 7 到 11 的平方幅度相加。 $f_0$ $N$ $N = \frac{4 F_s}{f_0}$ $N$ $4k+1 = 9$ $f0$

下面是一些示例代码，可以让您了解所讨论的算法的作用：

%% Generate an example signal
Fs = 5e3;                     % Sampling rate
f0 = 151;                     % Frequency of interest
t = (0:Fs-1)/Fs;              % Time samples
x = cos(2*pi*f0*t);           % Pure sinusoid @ f0
noise = 0.1*randn(size(x));   % Gaussian noise
y = (x>0)-(x<0) + noise;      % Introduce harmonics and noise

%% Compute energy in the harmonics of f0
N = round(4*Fs/f0);           % Determine record length
Y = fft(y(1:N))/N;            % Compute N-FFT (normalized)
Y2 = Y.*conj(Y);              % Magnitude squared
K = floor(N/2/4);             % Number of harmonic "bins"
for k = 1:K
    idx =  4*k+1+(-2:2);      % Five indices centered at 4k+1
    Pyk(k) = sum(Y2(idx));    % Sum of energies around 4k+1
end

%% Plot the results
[Pyy,Fyy] = pwelch(y,[],[],[],Fs);
plot(Fyy,10*log10(Pyy));
hold on;
F = linspace(0,Fs-Fs/N,N);
Fk = F(1+4*(1:K));
plot(Fk,10*log10(Pyk),'ro');
title('Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/rad/sample)')
data2 = sprintf('Energy in Harmonics of %d Hz',f0);
legend('Power Spectral Density Estimate', data2);
hold off;

来自 Matlab 代码的示例图

至于在数据中找到“签名”，您可以看到该算法可能对此非常有用。然而，非常重要的是要注意该算法假定了所讨论信号的一些基本频率。如果所讨论的信号具有固定周期，但该周期内信号的形状随时间而变化，则此算法将很有用。但是，如果频率随时间变化，则此算法不适合您的应用。

如果您只是想估计信号的频谱，那么此页面可能会有所帮助：FFT for Spectral Analysis (mathworks.com)

其它你可能感兴趣的问题

上一篇音乐噪音听起来像什么，是否有样本剪辑？下一篇关于医学成像的书籍