计算信号的 fft 频谱并获得正确的幅度

Question

我正在尝试一些相当简单的事情，即在具有固定长度的信号的频谱中获取正确的值。我不知道包含的频率，因此我无法相应地选择 fft 节点的数量。这给我留下了泄漏的问题。我试图通过使用 window 淡出我的时间信号的开始和结束来避免这种情况。因此，我在获得正确的振幅方面遇到了问题。这是我到目前为止所做的：

通过将频谱乘以一个因子
双边到单边频谱 -> 将频谱幅度乘以二
normalice 频谱振幅除以频率节点数'N'

只要 fft 节点的数量等于我的时间信号中的数量，它就可以正常工作。但是当 fft 节点的数量少于时间信号节点时，频谱的幅度就会出错。

我构建了一个计算三个光谱的小例子：

将 fft() 与 windowing
在没有 windowing
使用频谱图()

改变SampleTime可以改变时间信号的长度，改变N可以改变FFT节点数。使用以下设置，所有振幅都是正确的，但对于具有 Hanning window 的频谱。为什么差那么多？提前感谢您的帮助。

SampleTime=0.001;
Fs=1/SampleTime;
t=0:SampleTime:10;

y=2*sin(2*pi*10*t)+3*cos(2*pi*30*t); % Time Signal

fn = Fs/2; % Nyquistfrequency
N = 1001; % FFT-length (use N=2^x to use DFT-Algorithm)
df = Fs/N; % frequency resolution

%% Generating spectrum using fft()
win = hann(length(y))';  % generate Hanning window for Time Signal

ywin = y.*win*length(win)/sum(win); % amplitude correction to get correct values with window

Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window

amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/N);
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/N);
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

%% Generating spectrum using spectrogram()
noverlap = floor(0.5*N);
window = hann(N);
[B,F,T] = spectrogram(y,window,noverlap,N,Fs);
% Scaling
FFTScalingFactor = 2 * 1 / N; % 2 for Hanning window, 2^(-1/2) for RMS, 1 for peak
A = FFTScalingFactor * B;
B = FFTScalingFactor * abs(B);
B(2:end-1,:) = 2.*B(2:end-1,:); % single side spectrum with even block length

%% Plots

x_fn = 0 : df : fn-df; % Frequency vector

figure
subplot(2,2,1)
plot(x_fn,amplitudengangwin,'r')
hold on
plot(x_fn,amplitudengang)
legend('Hanning window','no window')
title('FFT()')
subplot(2,2,3)
plot(t,y)
title('Time signal')
subplot(2,2,2)
plot(F,B(:,1))
title('Spectrogram()')

Answer 1

你应该改变FFT的数量。

Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window

此处N为1001

这意味着当您对 windowed 信号执行 FFT 时，您只使用了所有数据样本中的前 1001 个样本（10001 个样本）。

您使用的 windowed 信号如下所示。

像下面这样更改您的代码

Hwin = fft(ywin, length(y)); % spectrum of windowed signal
H = fft(y,length(y));  % complex spectrum of signal without window
    amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/length(y));
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/length(y));
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

fftAll=[amplitudengang' amplitudengangwin'];

figure(1)
    plot(fftAll)
    xlim([0 500])

你会得到低于频谱的结果。

您定义了 N=1001，但您的数据长度为 10001。如果您只想使用 1001 个样本，则应将 window 长度设置为 1001。

如果您只想在数据中使用 1001 个样本，请按如下所示更改代码。

SampleTime=0.001;
Fs=1/SampleTime;
t=0:SampleTime:10;

y=2*sin(2*pi*10*t)+3*cos(2*pi*30*t); % Time Signal

fn = Fs/2; % Nyquistfrequency
N = 1001; % FFT-length (use N=2^x to use DFT-Algorithm)
df = Fs/N; % frequency resolution

%% Generating spectrum using fft()
win = hann(N)';  % generate Hanning window for Time Signal

figure(10)
plot(win)

ywin = y(1:N).*win*length(win)/sum(win); % amplitude correction to get correct values with window
figure(11)
plot(ywin)
xlim([0 1001])
   Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window
    amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/N);
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/N);
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

fftAll=[amplitudengang' amplitudengangwin'];

figure(1)
    plot(fftAll)
    xlim([0 50])

在此代码中，我将所有长度(y) 更改为 N。我像下面这样改变了 ywin。

ywin = y(1:N).*win*length(win)/sum(win);

结果将如下所示。

计算信号的 fft 频谱并获得正确的幅度

Calculating fft spectrum of a signal and obtain correct amplitudes

matlab

signal-processing

window

fft