Correct Amplitude for f=0 and fs/2

[chris-hld]

Spectrum is symmetrical around 0 and fs/2, therefore the corresponding amplitudes should not be scaled by factor 2. Try for f = 0: sig = ones(1,100); Try for f = fs/2: sig = repmat([0; 1], 22050, 1);

[fietew]

FYI: http://dsp.stackexchange.com/questions/19817/using-only-positive-frequencies-in-fourier-domain-how-will-it-affect-the-ifft

[fietew]

There are several bugs in the easyfft and easyifft even before your changes:

First question to @hagenw: What is the purpose of these functions? And what are the inputs and outputs of these function exactly (e.g. length/meaning w.r.t a "normal" DFT/IDFT)?

First example (even signal length):

SFS_start;
conf = SFS_config;

sig = [1 1 0 0];
[amp, pha, f] = easyfft(sig, conf);

sig_back = easyifft(amp, pha);

The "normal" DFT has frequency samples at f = [0, fs/4, fs/2, fs*3/4] or due to periodicity at f = [-fs/2 -fs/4 0 fs/4] . The spectrum values are unique at 0 and fs/2, while the values at -fs/4 and fs/4 have a complex conjugate symmetry. So the expected outputs of easyfft should have length 3 imho. However, the easyfft returns only 2 values, which is due to a wrong computation of the number of samples. It should be samples/2+1 or ceil((samples+1)/2) instead of ceil(samples/2).

Second example (odd signal length):

SFS_start;
conf = SFS_config;

sig = [1 0 0];
[amp, pha, f] = easyfft(sig, conf);

sig_back = easyifft(amp, pha);

The "normal" DFT has frequency samples at f = [0, fs/3, fs*2/3] or due to periodicity at f = [-fs/3 0 fs/3]. The spectrum values are unique only at 0. The outputs should have length 2 , which agrees with easyfft, but we could also change ceil(samples/2) to ceil((samples+1)/2) as both expression are equivalent for odd numbers.

The problem with the scaling of the amplitude is the following: If someone uses easyifft with a spectrum which he has generated by himself, he does not know, how he has to scale specific values of this spectrum in order to be valid for this function.

I would rather suggest to describe in the comments, that the functions do compute/accept the positive frequency part of the DFT spectrum. I think the scaling cannot be done in a sensible manner.

[hagenw]

Thanks for the review @fietew. At the moment I use easyfft() as a convinient way to create a FFT plot as it returns the positive spectrum and the corresponding frequency axis. When I wrote the function (I guess 10 years ago) I used it also to do some manipulation of signals in their amplitude-spectrum.

[fietew]

So is a scaling even necessary?

[chris-hld]

Yes, this was the second thing I wanted to adress. For the case where the signal length is even, the spectrum must contain [0, fs/2], which is does now. Whereas the ifft gets [0, fs[, after the the spectrum mirroring, due to the 2pi periodicity. This should work now.

I also added the assumption for the ifft, that the input is complex conjugate, in case there are round-off errors which will cause that the output are not only real values. Unfortunately I just figured out, that octave does not provide this option. So maybe we stick withreal(ifft(spec)).

I'm not sure how to handle odd number signals. However, I agree that in your example the output should have length 2. Which causes the ifft to generate a signal length 4 - with amplitude sig_back = [1.33 0 0 0] instead of sig = [1 0 0], because of the wrong amount of samples. But as mentioned in your link above, this error is very low for larger signals. So should we somehow provide a warning instead?

[fietew]

Whereas the ifft gets [0, fs[, after the the spectrum mirroring, due to the 2pi periodicity.

No, that is not correct in general. In easyifft two case have to be considered:

1. If the number of frequency bins is even, e.g. at f = [0 f/3], then the DFT stems from a signal of odd length. The correct mirroring is then done by `amp = [amp, amp(end:-1:2)].
2. If the number of frequency bins is odd, e.g. at f = [0 f/4 f/2], then the DFT stems from a signal of even length. The correct mirroring is then done by `amp = [amp, amp(end-1:-1:2)].

IDFT(DFT(sig)) must have the same length as sig

[chris-hld]

True that, thanks! Maybe its best to treat each case individually. Its easier to understand and we can provide the correct scaling.

[fietew]

I'm still not convinced that this scaling, in general, is a good idea. Imho, it has no meaning/reason in the sense of signal processing.

[fs446]

Scaling gives you the correct amplitudes of DC (pi), the first DFT harmonic and the DFT harmonic at fs/2 (amplitude 1) in this example:

N = 8; x = pi + sin(2_pi/N1[0:N-1])+cos(2_pi/NN/2[0:N-1]); X= fft(x); Xabs = 2/N_abs(X); Xabs(1) = 1/2_Xabs(1); if mod(N,2)==0 Xabs(N/2+1) = 1/2*Xabs(N/2+1); end

This is required when you want to have correct plots or you want to synthesize signals from DFT spectra, where DC and fs/2 components are important.

Best Frank

On 16 Sep 2016, at 2:49 PM, Fiete Winter notifications@github.com wrote:

I'm still not convinced that this scaling, in general, is a good idea. Imho, it has no meaning/reason in the sense of signal processing.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[chris-hld]

Yep, thanks Frank! But something doesn't work for odd signals: An odd signal with length 9, returns 5 frequency bins [0 4/9]*fs, which is correct. But this causes the ifft to return 8 values... What am I missing?!

Maybe the right way is testing whether the frequency vector contains fs/2 or not?

[fietew]

Oh, I think my statement about the relation between signal lengths and bin numbers was wrong. You can't separate e.g. f = [0 f/5 f*2/5] and f = [0 f/4 f/2]

[chris-hld]

As long as the sample rate is not incredibly high, the last entry of the frequency vector can be used to decide for the length of the IFFT output. Therefore both the f vector and conf.fs is needed as function arguments. Any smarter ideas for this? The function passed all my simple tests now.

[hagenw]

If no one else has any complains with the current status of this proposed changes I will merge them at the end of this week.

[fietew]

A user, who is not aware, that he has to multiply the value of the spectrum at 0 and fs/2 by 0.5 (or the other values by 2), will get a incorrect result. Is that correct?

[chris-hld]

No, thats exactly the point! This is unnecessary now, and amplitudes are consistent for both time and frequency domain! Have a look also at the validation example. It's obvious there. I added an example with fs/2 as well.

[fietew]

I'm aware, that a spectrum created with easyfft is compatible with easyifft and easyifft(easyfft(x)) yields x again. This is out of question.

But what if a user defines the Fourier spectrum manually for the positive frequency axis and wants to get the correct (real) signal. Normally, he would use the ifft( , 'symmetric') from MATLAB for that, but he has to be careful, as he has to append the right number of zeros to the positive spectrum. So now he thinks that easyifft() eases things up for him, but he can't use this function straightforwardly.

conf = SFS_config;
fs = conf.fs;

% desired spectrum for positive frequency axis
X = [1, 1, 1, 1, 1];
f = linspace(0,fs/2, 5);
% conventional way
Xfft = [1, 1, 1, 1, 1, 0, 0, 0];
xfft = ifft(Xfft, 'symmetric');
% way with easyfft?
xeasy = easyifft(abs(X), angle(X), f, conf);

The signals xfft and xeasy are not the same. Again the question arises: For what purpose do we use the functions easyfft and easyifft? The description easyifft says "EASYIFFT calculates the inverse FFT", but it does that only if the user is aware of the type of spectrum he has to provide.

[fietew]

see https://en.wikipedia.org/wiki/Fourier_series#Definition for explanation and nomenclacture

fft(x) and easyfft(x) essentially provide equivalent coefficients of a Fourier series of the discrete, real input signal x. Note, that the wiki entry assumes a continuous signal x so that the coefficients are not periodic w.r.t to n. While fft provides the complex coefficients cn=c*-n for the complex notation of the Fourier series, easyfft provides A_n and phi_n (n>=0) for the amplitude-phase-notation of it.

Imho, it would be way more intuitive to provide the c_n (n>=0), i.e. the coefficients of fft for the positive frequencies than providing A_n and phi_n especially because the function is called easyFFT.

[chris-hld]

I agree on updating the docstring of easyifft. Apparently, the option 'symmetric' does more than documented. The docu says it forces the output to be real, when your spectrum has round off errors/ is not conjugate complex anymore.

This option is useful when X is not exactly conjugate symmetric, merely because of round-off error

Actually, easyifft should just be the complement to easyfft. I mainly know cases, where you are interested in absolute amplitude and phase behaviour of a system, maybe manipulate that, and then create a time signal again. That is, as you said, the notation

A_n and phi_n (n>=0)

Otherwise, if there is already a complex spectrum, there is no need for easyifft imho. But I'm interested, why your version creates a perfect dirac vs. easyifft creates a dirac with leakage, that's not obviuos to me right now...

[hagenw]

I would not be in favor of getting back c_n, as the main reason of creating easyfft() was to get A_n and phi_n without doing anything. We can also consider to rename it, or to remove easyifft() if you think that people will use it in a wrong way?

[fietew]

@chris-hld: As you already said, the behaviour of 'symmetric' is not documented very well. Internally, the function throws the negative frequency samples away and computes a real-valued IFFT from the positive frequency values, which is a little bit faster than the normal FFT and does not provide small imaginary parts. However, the user has to provide to correct length of the whole spectrum, because the algorithm has to decide if it has to keep the sample for fs/2. As the negative frequency are thrown away, the user could append any values.

easyifft assumes that all frequency samples except of f=0 and f=fs/2 are weigthed by two and that the normalisation to the number of samples is already done by easyfft. This should give the correct result:

conf = SFS_config;
fs = conf.fs;

% desired spectrum for positive frequency axis
X = [1, 1, 1, 1, 1];
f = linspace(0,fs/2, 5);
% conventional way
Xfft = [1, 1, 1, 1, 1, 0, 0, 0];
xfft = ifft(Xfft, 'symmetric');
% way with easyfft?
X = [1, 2, 2, 2, 1]/8;
xeasy = easyifft(abs(X), angle(X), f, conf);

@hagenw: Personally, I use easyfft() only for quick plotting. I agree, for this case the changes in this pull request make sense. Imho, we should keep easyifft(), but however clearly state, that both functions are not compatible with a "cross-use" together with fft() and ifft(). I would prefer to rename both functions in order make the difference ever clearer. I also had a quick chat with @spors and he agrees, that the output of easyfft is not a Fourier transform anymore.

[narahahn]

Hi guys, what an interesting discussion. I'd like to share my thoughts on this issue.

Since easyfft and easyifft have fft in their names, one may expect that c_n = A_n * exp(1i*phase) for n>0. But this is not the case. In DFT, the amplitude response at 0 and fs/2 are scaled. This is not intuitive but this is how DFT is. In my opinion, it would be more clear to maintain the c_n = A_n * exp(1i*phase) relation. If we want to keep these functions as they are, it would be better to give new names, like 'freqresp`.

easyfft and easyifft take fs as an input argument but normally we do not need fs for DFT. It come into play if we want to plot a frequency response. @hagenw : are these functions mainly (if not only) for plotting frequency responses? Imho, if so, wouldn't it be better to use an extra function for plotting frequency responses and consider the scaling there? Just my opinion:-)

If a manipulation of an amplitude spectrum is required, we can do it for c_n as well. If we do so, we can also apply phase response simultaneously, instead of multiplying a scale to amplitude and adding a phase to phase. Scaling at 0 and fs/2 is not an issue here. If we want an boost of let say 3dB at 0 or fs/2 it will just work. A phase shift at those frequency bins will be automatically ignored after real(ifft(X)) or `ifft(X,'symmetric').

[hagenw]

Hi @narahahn, thanks for your input. I personally used the function for these two purposes:

1. Easy way to plot the frequency response
2. Grab a signal, transform it with easyfft, manipulate its amplitude spectrum, and transform it back with easyifft (I have to admit, that I haven't done this for a long time, but during my psychoacoustics days in Oldenburg, this was done quite often).

[narahahn]

Sorry I was wrong. We have to be careful with the phase response at 0 and fs/2. It will not be ignored.

[chris-hld]

Alright guys, thanks for all the great input! How should we proceed? I agree, the naming easyfft is not very suitable. It is not a "pure" fourier transform anymore, since it is extended quite a bit towards describing the spectrum. So I propose the following:

• Rename easyfft to something like showspec (really not convinced with that yet).
• Optimize for instant plots of amplitude and phase spectrum
• Rename easyifft and make clear, it's only the counterpart for getting back a time signal

[hagenw]

The name get_spectrum is fine with me, but retrieve_signal is not ideal. The biggest problem in my opinion is the now missing direct link to the other function. Unfortunately, I have no better proposal at the moment, but I will try to find one.

[chris-hld]

BTW, the previous discussion probably makes retrieve_signal appear less generic than it actually is. In usual use cases it performs just as expected. The case of creating a dirac is quite specific, since the spectrum must contain every (two-sided) frequency with the same amplitude (summing up to 1), although only a single-sided spectrum is given. Therefore, the energy of discarded frequency bins must be contained in the mirrored bins (that's the factor 2). For creating a 500 Hz cosine from a single-sided spectrum, try:

SFS_start; conf = SFS_config
f = [0:1:22050];
ampl = zeros(1,22051);
phase = zeros(1,22051);
ampl(501) = 1;
signal = retrieve_signal(ampl, phase, f, conf);
get_spectrum(signal, conf)

Which is exactly what I would expect.