フィルターのクリア

Why does frequency spectrum gets flipped in STFT compared to the one generated from FFT?

16 ビュー (過去 30 日間)
Edward
Edward 2023 年 8 月 30 日
編集済み: Edward 2023 年 9 月 1 日
I picked out a signal of 200 temporal samples and processed it separately with STFT and FFT. The data is attached as .mat file.
The code for STFT is listed below (the first dimension of Matrix has 200 temporal samples and I plot the spectrogram of first signal segment):
[s,f,t] = stft(Matrix',5000/3,Window=rectwin(32),OverlapLength=ceil(0.8*32),FFTLength=200); % Full map
s_amp = abs(s);
figure;plot(squeeze(s_amp(:,1,1)));
The code for FFT is shown here (I get the first 32 temporal samples and do the FFT, then the spectrum is zero centered):
FreqSpec_temp = fft(squeeze(Matrix(1,1:32)),[]);
FreqSpec_shift = fftshift(FreqSpec_temp);
figure;plot(abs(FreqSpec_shift));
The plots of STFT (left) and FFT (right) are displayed. One can easily spot out that the frequency spectrum of STFT is horizontally flipped compared to the spectrum of FFT. Why is this the case? Which method should I use to get an accurate frequency spectrum? I will need to calculate the mean frequency of the signal based on the spectrum.
Besides, how could the spectrum from STFT has 200 points? How does it work out the interpolation? Thank you for your help!
  2 件のコメント
Paul
Paul 2023 年 8 月 31 日
編集済み: Paul 2023 年 8 月 31 日
Hi Edward,
It will be easier for someone to help if you save M in a .mat file and upload it using the paperclip icon in the Insert menu so that the code and figures can be recreated.
"how could the spectrum from STFT has 200 points? "
Edward
Edward 2023 年 8 月 31 日
Hi Paul,
Thank you for your comments. I just uploaded the .mat file for you to try out.
For FFT, when it takes in 32 temporal samples, the frequency spectrum will also have 32 points for representing the spectrum. But for STFT, it takes in 32 temporal samples for processing 1 temporal segment, while the frequency spectrum has 200 points. What I wonder was certain interpolation was done in generating the frequency spectrum, but the process was not specified in the page of STFT. I would like to know what kind of interpolation was used to get 200 points from an input of 32 points.

サインインしてコメントする。

採用された回答

Paul
Paul 2023 年 8 月 31 日
Using random data
rng(100);
Matrix = rand(1000,1);
Compute the stft
[s,f,t] = stft(Matrix,5000/3,Window=rectwin(32),OverlapLength=ceil(0.8*32),FFTLength=200); % Full map
s_amp = abs(s);
plot amplitude vs frequency
figure;plot(f,squeeze(s_amp(:,1,1,1)));
Compute the fft, zero pad to 200 to match FFTLength
FreqSpec_temp = fft(squeeze(Matrix(1:32,1,1)),200);
FreqSpec_shift = fftshift(FreqSpec_temp);
Plot amplitude vs frequency, matches first column of stft
hold on
plot((-100:99)/200*5000/3,abs(FreqSpec_shift),'o');
  3 件のコメント
Paul
Paul 2023 年 8 月 31 日
編集済み: Paul 2023 年 9 月 1 日
I think the issue may have been that the code was using the conjugate transpse (ctranspose, ') on the input to stft, but was not taking the conjugate on the input to fft. I've corrected that below, but it might not be what you want. I suspect that you really want to use transpose (transpose, .') on input to stft (and not take the conjugate on input to fft).
Because the input to stft specifies FFTLength=200, the stft will have 200 points covering the "centered" frequency interval, which is the default and used here. I'd refer you to the doc page, but, frankly, the doc page seems to be in error in defining what that frequency interval is. Also, stft with the "centered" option uses a different convention than fftshift, which is why I didn't use the f output from stft for plotting the fft.
In short, using a 200-point fft on a 32-point input sequence yields 200 frequency domain samples of the Discrete Time Fourier Transform (DTFT) of the 32-point input sequence, so it's not really interpolating, per se (though it's typically called "frequency interpolation"). I'm sure that readily available references can be found that discuss the relationship between the DTFT and the Discrete Fourier Transform (DFT) that is computed with fft.
load('data.mat','Matrix')
whos
Name Size Bytes Class Attributes Matrix 2x200 6400 double complex ans 1x33 66 char cmdout 1x33 66 char
This line is copy/pasted from the Question. Note the use of ctranspose, ' on the first argument, which is the conjugate transpose, which is probably not what's desired.
[s,f,t] = stft(Matrix',5000/3,Window=rectwin(32),OverlapLength=ceil(0.8*32),FFTLength=200); % Full map
s_amp = abs(s);
figure;plot(f,squeeze(s_amp(:,1,1)));
Because the input Matrix to stft was conjugated, we need to do the same on the input to fft. And need to zero-pad the fft to 200 points to match the FFTLength used in stft
FreqSpec_temp = fft(squeeze(conj(Matrix(1,1:32))),200);
FreqSpec_shift = fftshift(FreqSpec_temp);
Plot amplitude vs frequency, matches first column of stft
hold on
plot((-100:99)/200*5000/3,abs(FreqSpec_shift),'o');
If we don't zero pad the input to fft, the 32 output samples cover the same frequency range, but have larger spacing
FreqSpec_temp = fft(squeeze(conj(Matrix(1,1:32))),[]);
FreqSpec_shift = fftshift(FreqSpec_temp);
figure
plot(f,squeeze(s_amp(:,1,1)));
hold on
plot((-16:15)/32*5000/3,abs(FreqSpec_shift),'o');
Note that the frequency points for the two curves aren't the same, which is why the there is an apparent difference around f = -100. But that's just because of where the respective frequency samples lie
figure
plot(f,squeeze(s_amp(:,1,1)),'-o');
hold on
plot((-16:15)/32*5000/3,abs(FreqSpec_shift),'-o');
xlim([-200 0])
Edward
Edward 2023 年 9 月 1 日
編集済み: Edward 2023 年 9 月 1 日
Thank you, Paul! Yes, you're right. I wanted to use transpose (.') only, but I didn't realize ' actually represents as conjugate transpose. It was my coding mistake. And thank you for explaining the difference between the 2 curves!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeDescriptive Statistics についてさらに検索

製品


リリース

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by