이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
Correct choice of one sided frequency axis after fft
조회 수: 8 (최근 30일)
이전 댓글 표시
There are a lot of queries on fft frequency. I guess the following point not discussed anywhere explicitly (at least at one place). Hope someone can provide an insight here.
If we have and even number of data points, N=10, the fft output arranges the data as
fft = [c0, c1, c2, c3, c4, c-5, c-4, c-3, c-2, c-1], where the numbers are subscripts corresponding to positive and negative frequencies. I read somewhere that MATLAB calculates the negative coefficient first, hence we have c-5 but not c5. The author did not explain the reason.
Point no. 1, that the c values are not symmetric.
When we wish to make two-sided frequency spectrum, the frequency axis ranges from [-(N/2): (N/2)-1]*Fs/ N. Fs is the sampling rate, N is the number of even data points.
If we wish to make a one-sided positive frequency spectrum, should we choose
A) [0:(N/2)]*Fs/N and ignore the fact the we are using the values corresponding to the negative frequency axis, given that the data is a real number and it is just a mirror image.
B) [0: (N/2)-1]*Fs/N represents the true positive frequency axis?
If Fs= 250 Hz, the true positive frequency axis will end at 124.9980 Hz
If we happen to choose the negative frequency axis values and ignore the frequency sign, the frequency axis ends at 125 Hz exactly.
The same data when plotted in Origin ends the frequency axis at 125 Hz when plotted single sidedly.
Which approach is rigorously correct? Thanks.
댓글 수: 21
dpb
2019년 7월 7일
Why would you choose anything but the actual DC to max positive frequency that is actually calculated? That there is one sample short of Fs/2 is just a fact of "the way things are".
The time history of N points beginning at T=0 is one element short of N*dt as well...maybe that realization will help to understand the "why" of it all...
FW
2019년 7월 7일
Why would you choose anything but the actual DC to max positive frequency that is actually calculated? That there is one sample short of Fs/2 is just a fact of "the way things are".
Yes it is. The key question is a conceptual problem rather. I am not worrying about the positive frequency axis ending at 124.9980 Hz rather than 125 Hz. but why does OriginPro have positive the frequency axis ending at 125 Hz exactly and MATLAB doesn't?
Are they including fft = [c0, c-1, c-2, c-3, c-4, c-5]?
dpb
2019년 7월 7일
I have no idea; you'll have to ask them or see what the values are that they have used are empirically and determine.
This is the wrong forum for Origin-specific queries; this is Matlab.
FW
2019년 7월 7일
Is it officially documented somewhere by Matlab that the fft out is in this order as shown below? This problem is very much relevant to Matlab's way of indexing, Origin was just used for comparison.
fft = [c0, c1, c2, c3, c4, c-5, c-4, c-3, c-2, c-1]
Thanks.
Bruno Luong
2019년 7월 7일
Discrete Fourier Transform of Vector
Y = fft(X) and X = ifft(Y) implement the Fourier transform and inverse Fourier transform, respectively. For X and Y of length n, these transforms are defined as follows:
Y(k)=SUM j=1...n X(j) W_n(j−1)(k−1)
X(j)=1/n SUM k=1...n Y(k) W{n}-(j−1)(k−1),
where
W_n=e(−2πi)/n
is one of n roots of unity.
The problem often arises when people try to interpret the formula as frequency and they often mess up or confuses themselve.
FW
2019년 7월 7일
Bruno, I am talking about the indexing and ordering of the fft output by Matlab. I have seen this Matlab page you linked. Also see a similar discussion here, and the problem is perhaps not as trivial.
Bruno Luong
2019년 7월 7일
The index location / ordering is fully defined in this formula. Don't you see it ?
What is not trivial for you is when you talking about c5 and c-5 as if they are two different things, they actually represents the same function if you look at for formula and it located at y(6) (and it's Nyquis frequency). Period.
dpb
2019년 7월 8일
As the doc says, ML uses the fftw implementation so there's nothing unique about it other than ML does store all arrays with one-based indexing. That's simply bookkeeping.
ML returns the DC component in the first array position, then the fix(N/2) positive followed by the same number of negative frequencies in obverse order. If N odd, then the returned spectrum is symmetric about DC.
When N even, it is simply not possible to have a DC component and N/2 positive and negative frequency components...just the way it is. fftw will have the same "problem" and Origin will too unless they do something like either augment or truncate the input vector to be odd internally or simply add another element or somesuch on the output.
Bruno Luong
2019년 7월 8일
편집: Bruno Luong
2019년 7월 8일
The DFT formula never tells which Y elements corresponds to "positive" or "negative" frequency. Often user interprets them as such but DFT never state signed frequency, since one can add an (k'*n), where k' is an arbitrary integer in "frequency" to give the same result.
In the example of n=10, the c-5, c5 are the same coefficient of the fast oscillated sequence (1,-1,1,-1, ....-1) and I can arbitrary view it as c25 or c1005, nothing is wrong with it since it gives the same result.
dpb
2019년 7월 8일
편집: dpb
2019년 7월 8일
True, but if one is looking at it from the viewpoint of baseband frequency analysis as is the OP, the ML storage order begins with the DC component up, then the obverse.
>> x=randn(1,11);
>> y=fft(x.')
y =
1.9397 + 0.0000i
2.7291 + 0.7234i
-0.6829 + 0.2768i
4.2973 + 1.0150i
3.0071 + 2.0347i
-0.0981 + 2.2543i
-0.0981 - 2.2543i
3.0071 - 2.0347i
4.2973 - 1.0150i
-0.6829 - 0.2768i
2.7291 - 0.7234i
>> mean(x)
ans =
0.1763
>> y/11
ans =
0.1763 + 0.0000i
0.2481 + 0.0658i
-0.0621 + 0.0252i
...
The first element returned by FFT after normalization is the mean of the input--the DC component if doing typical signal processing of a measured input.
FW
2019년 7월 8일
Hello dbp, I did a quick check with OriginPro double sided spectrum. The frequency scale is automatic there when sampling rate is provided. At 250 Hz, with OriginPro, the negative freq. ends at -124.9980 Hz, and positive frequency ends at +125 Hz. However, if we follow the Matlab output order, the positive frequency ends at +124.9980 Hz and negative frequency end is -125 Hz. It seems like a matter of convention and both approaches by Matlab or OriginPro are "mathematically" correct.
This was my Matlab code
t=[0:1/250:10.004]'; % ensuring even number output
Fs= 250; % Hz units
N= length(t);
g=normpdf(t, 5, 0.2); % Gaussian peak of unit area
f_Hz=[(-N/2):(N/2)-1]'*Fs/N;
G_double=fft(g);
fftshift_G = fftshift(G_double); % fftshift rearranges the vectors so that zero frequency is in the center
Y=abs(fftshift_G);
plot(f_Hz, Y, '.')
Bruno Luong
2019년 7월 8일
편집: Bruno Luong
2019년 7월 8일
"if we follow the Matlab output order, the positive frequency ends at +124.9980 Hz and negative frequency end is -125 Hz."
Again no : MATLAB FFT doc/help doesn't never ever tell such thing. This is an interpretation by you (or whoever you read and believe). We can select either convention of frequency intervals
[-124.9980,125] Hz, or
[-125,124.998] Hz.
Both are equally correct, or even an odd choice of frequency interval:
[-125,124.998] + 250 Hz.
Odd but still correct.
FW
2019년 7월 8일
Again no : MATLAB FFT doc/help doesn't never ever tell such thing. This is an interpretation by you (or whoever you read and believe). We can select either convention of frequency intervals
[-124.9980,125] Hz, or
[-125,124.998] Hz.
This makes much more sense. My source of information was a primer on FFT for Physicists. The author does not cite anything. Basically you are saying that there is nothing in MATLAB which supports form of the table. This is the table I was following (full pdf attached), see page 10.
Bruno Luong
2019년 7월 8일
The paper cites about the choice of FFTSHIFT applied on even-data spectrum.
This choice is arbitrary. User can take it or leave it, in the later case placing the nyquise as positive if he/she prefers (that makes compliance with OriginPro).
dpb
2019년 7월 8일
편집: dpb
2019년 7월 9일
I think your author goes out of his way to confound what is really quite simple... :)
I'd basically ignore the table as being of no benefit for any purpose under the sun.
As noted earlier, ML returns DC component in first array location; whether there are even or odd number in the series determines whether the other components are or are not symmetric in pairs or whether there's an "odd man out" when even.
As Bruno says, depending on your purpose you can either fftshift or not depending upon objective at hand. I rarely bother, simply take the first half as I don't recall in 40 years a need for anything but the one-sided spectrum. Others' needs clearly may differ... :)
FW
2019년 7월 9일
Thanks dpb, I guess it all boils down to convention. This FFT primer was pretty good and I thought this is the way Matlab arranges the output (as shown in the table). However I could not find any other page which confirmed it. It is now clear that this the author's own interpretation, which is not wrong either. I am not a mathematican nor a signal processing engineer. As a chemist, I was trying to implement a fft step in resolution enhancement protocols for a spectroscopic signals. The take home message is that fft output, say N=6, it is
fft = [c0, c1, c2, c3, c4, c5], where the number k refers to N-1 values.
for which we can construct the frequency axis as
-Fs/2: (Fs/2-Fs/N)
or -(Fs/2-Fs/N): Fs/2 is completely up to us, and both frequency axes are correct.
dpb
2019년 7월 9일
>> n=0:9;
>> [n;fftshift(n)]
ans =
0 1 2 3 4 5 6 7 8 9
5 6 7 8 9 0 1 2 3 4
>> n=n+1;
>> [n;fftshift(n);fftshift(n)-1]
ans =
1 2 3 4 5 6 7 8 9 10
6 7 8 9 10 1 2 3 4 5
5 6 7 8 9 0 1 2 3 4
>> n=1:9;
>> [n;fftshift(n);fftshift(n)-1]
ans =
1 2 3 4 5 6 7 8 9
6 7 8 9 1 2 3 4 5
5 6 7 8 0 1 2 3 4
>>
Bruno Luong
2019년 7월 9일
편집: Bruno Luong
2019년 7월 9일
@dpb: Don't understand what you are at by doing
fftshift(n)-1
???
Don't you want to show something along this line:
ffts_l=@(x) circshift(x,floor((size(x)-1)/2)) % OriginPro compatible
ffts_r=@(x) circshift(x,floor((size(x))/2)); % equilalent to MATLAB fftshift
>> x=0:8;
>> [x; ffts_l(x); ffts_r(x)]
ans =
0 1 2 3 4 5 6 7 8
5 6 7 8 0 1 2 3 4
5 6 7 8 0 1 2 3 4
>> x=0:9;
>> [x; ffts_l(x); ffts_r(x)]
ans =
0 1 2 3 4 5 6 7 8 9
6 7 8 9 0 1 2 3 4 5
5 6 7 8 9 0 1 2 3 4
dpb
2019년 7월 9일
편집: dpb
2019년 7월 9일
Bruno, the n-1 is just the offset difference between 1- and 0-based arrays--a trivial difference but seems to throw many into complete tizzy.
Your latter illustration is also pertinent, I was just too lazy to write the extra code figuring simply showing OP the indexing is the root answer of how Matlab returns the order and one can build the quoted author's table from it by inspection.
Bruno Luong
2019년 7월 9일
편집: Bruno Luong
2019년 7월 9일
"difference between 1- and 0-based arrays"
It escapes me about who was discussing of 1/0 base indexing in this thread.
As I understood the discussion is around the FFT output coefficients c_i and its subindex i which makes the correspondance to the frequency:
f(i) := Fs * i / n.
and whereas i should vary beween
[-floor((n-1)/2) : ceil((n-1)/2)] % OriginPro
or
[-ceil((n-1)/2) : floor((n-1)/2)] % MATLAB FFTSHIFT
(As we discuss, both are mathematically correct and it's a matter of interpretation or arbitrary choice)
Bjorn Gustavsson
2019년 7월 9일
@dpb (on 7 Jul 2019 at 14:50): One "common" use of "not choosing" base-band is that one has undersampled a high-frequency signal with signal at frequencies inside the N-th Nyqvist zone. The spectra of that signal can then be accurately determined with sampling-rates much lower than the centre-frequency of the band.
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Spectral Measurements에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)