What is the proper way to compute a real-valued time series given a continuous $1/\sqrt{\omega}$ spectrum?

Question

I have never fully been able to wrap my head around Fourier transforms, so I apologize if what I am trying to do is trivial or violates basic theory in some way.

What I have is a "made up" frequency spectrum (well, it comes from somewhere but not relevant for my problem). The spectrum has the form:

$$|G(\omega)| = \frac{a}{\sqrt{b\omega }}$$

where $a$ and $b$ are arbitrary (real) constants for my particular application. Note that I only show the spectrum but the real complex-valued Fourier spectrum would be: $$G(\omega) = \frac{a}{\sqrt{b\omega }}e^{i\phi(\omega )}$$

For my purposes, $\phi(\omega)$ can be anything and is not physically constrained. As a starting point I've just set $\phi(\omega)=0$ so that the spectrum is real-valued. Sampling rate (or frequency spacing) can be arbitrary as can the length of signal.

What I want to do is find a time series, $g(t)$ that is the inverse Fourier transform of $G(\omega)$.

What's the easiest way to do this?

Thus far, I've done it an exceedingly "dumb" way that probably makes many specialists in time series processing cringe. I've basically just digitized the frequency spectrum like a comb filter and then calculated the sinusoids at each spike in the comb and then added them together. Basically treating each individual digitized frequency as a delta spike (implemented in MATLAB below), using the most basic defintion of the Fourier transform:

$$ g(t) = \int_{-\infty}^{\infty} A(\omega)e^{i2\pi\omega t} d\omega$$ The resulting time series I get actually kind of looks "okay" on first glance. However, if I try to perform a standard MATLAB fft of this time series, it returns complete garbage (see figures below)

Code:

N = 1000; %number of discrete frequencies to digitize the continuous spectrum
a = 1.2566; %constant
b = 0.00078956; %constant

w = 10.^linspace(0,-4,N); %get logarithmically spaced frequencies omega

%The spectrum
G = a./sqrt(b*w);


M = 100000; %length of time series
fs = 1; %assuming a 1 Hz sample rate

t = linspace(0,fs*M,M)';

%the time series is computed as a "dumb" Fourier transform which assumes
%that each discrete frequency is a delta spike respresenting a time-domain
%sinusoid. Just sum it all together:
g = sum(G.*exp(1i*2*pi*w.*t),2); 


%Now perform fft on the time series you just made
pad = 10000;
xt = cat(1,ones(pad,1).*g(1),g,ones(pad,1).*g(end)); %add some padding cells to the time series
X = fft(xt); 

%Calculate frequency bins
df = fs/length(xt);
fAxis = (0:df:(fs-df)) - (fs-mod(M,2)*df)/2;

%Plotting
subplot(2,1,1) %Plot spectrum
loglog(w,G,'*r'); hold on; grid on
loglog(fAxis,abs(X),'-k'); 
xlabel('Frequency (Hz)');
ylabel('G(\omega)')
set(gca,'FontSize',14)
title('Spectrum G(\omega)')
legend('"True" Spectrum','FFT of Time Series')

subplot(2,1,2); %Plot time series
plot(t/3600,real(g)); hold on;
xlabel('Time (hours)')
ylabel('g(t)'); grid on
set(gca,'FontSize',14)
title('Inverse "Fourier Transform" of Spectrum')

Result:

Any help is appreciated!

Answer 1

Bottom Line

• OP's approach is closest to the "Frequency Sampling Method"
• Use a linear frequency grid, not logarithmic.
• Populate the spectrum on both the positive an negative frequency axis, with complex conjugate symmetry for real time samples.
• (As MattL pointed out): bound the DC bin to a finite value.

The OP's used logarithmically spaced positive axis only frequencies but proceeded with what is equivalent to a Discrete Fourier Series reconstruction where linearly spaced frequencies should be used. Using linear spaced frequencies that include the positive and negative frequencies, and using the inverse FFT directly to estimate the time domain response results in the "Frequency Sampling" method. The plot below shows the achievable result of the recreated spectrum using 1001 samples as the OP had done. Increasing the number of samples extends the accuracy to even lower frequencies.

The remaining details in this post should clear up any disconnects with the OP's initial intuition with the approach that was attempted.

MORE DETAILS

The main points of this post, including what was wrong with the OP's initial attempt is concluded above. The remaining part of this post details the deeper considerations needed to do the Frequency Sampling Method properly for any generic spectrum (real or complex), and then as an interesting study I included a direct comparison to the result for the Windowing Method that MattL derived in his answer (initially I expected the Frequency Sampling solution using the same number of samples to be worst without increasing the density of samples in frequency and wanted to demonstrate that, but that didn't turn out to be the case). That later section also summarizes the differences and approach of the two methods.

The Frequency Sampling method is simple and quite effective if you are able to sample the frequency over the spectrum of interest with fine enough spacing to avoid time domain aliasing distortion. To proceed, sample the desired spectrum using both its magnitude and phase, and then use an inverse FFT calculation to determine the time domain samples. The finer the spacing in frequency that is used, the longer the resulting time domain response will be, which will further reduce deleterious time-domain aliasing effects. This resulting time response can be truncated to a smaller duration, and then this result windowed to reduce the time-domain truncation effects. Without such windowing the time-domain truncation will result in ripples in the frequency response. Extending the sampling rate higher (most simply by zero-padding the target spectrum) is also an effective approach to interpolate more samples in the resulting time domain response if needed.

There are some details that are important to be aware of:

Detail 1 With using the inverse FFT directly, the sampled frequency domain must extend from DC to 1 sample less than what will be the sampling rate for the resulting time domain signal. This is to properly fill the values that essentially are the FFT result: $N$ samples in time are $N$ samples in frequency, the first FFT bin with $k=0$ is at "DC" and the last FFT bin with $k=N-1$ is at the frequency that is nearly the sampling rate (if there was one more bin at $k=N$, it would be at the sampling rate). The upper half of this spectrum is equally the "negative frequencies" due to periodicity of the FFT. Assuming a real time domain response, each negative frequency will be the complex conjugate of each positive frequency.

Detail 2 Using an odd number of samples simplifies the process since we can fill the FFT directly with $k=0$ as the "DC" bin, and the next $(N-1)/2$ samples as the "positive frequencies", and then the next $(N-1)/2$ samples after that as the "negative frequencies" and symmetry in the positive and negative frequency is maintained. (Using an even number of samples is done with further treatment of the "Nyquist bin" that is at $N/2$ in this case, which I won't detail here but that processing is covered in other posts [link here would be helpful]). As a simple example, the spectrum given by samples [1, 2, 3, 3, 2] would be considered symmetric and result in an all real inverse FFT: Here the DC bin has the value 1 associated with $k=0$, the positive frequencies are given as [2, 3] associated with $k=1,2$, and the negative frequencies are given as [3, 2] associated with $k=3,4$, where the two values given as 2 correspond to the two closest frequencies to the DC bin.

Detail 3 As MattL pointed out in the comments, the OP's spectrum will go to infinity at DC. The solution is to bound the DC bin to the value for the bin at $k=1$ (which for the real spectrum is also the bin at $k=N-1$). For this particular spectrum, I did not see a lot of sensitivity to the resulting spectrum by modifying the DC bin over a fairly wide range, but reducing it's value is intuitive given the effective first order high pass that is given by finite observation times (if a more detailed explanation of this is needed, that is worthy of a completely separate question).

Detail 4 A zero phase spectrum will result in a non-causal time domain response which is resolved by circularly rotating the result to the right over half its length (using ifftshift on the non-causal inverse FFT result will do this directly). That would be the simplest approach to take, but this can also be done by adding a negative linear phase in the frequency domain (which is a delay in time) prior to taking the inverse FFT, sufficient to circularly shift an odd length inverse FFT to the right $(N+1)/2$ samples. Also for cases where complex conjugate symmetry is not maintained in the frequency domain between the first half of the spectrum and the second half (between the positive and negative frequencies), the resulting time domain signal will be complex. (In some applications we want this).

Detail 5 When using a real and symmetric spectrum, the time domain response itself will be real and symmetric. If truncating this time response to shorten the filter length (and windowing), the truncation should be done to keep the center of the response to maintain the time domain symmetry. This means remove the same number of samples from the beginning and the end of the time domain response, and then window that result.

If those details are properly maintained when going from frequency to time using the inverse FFT, then the resulting time domain response will have a continuous frequency response (The Discrete-Time Fourier Transform or DTFT) with close match to the original desired spectrum.

In summary the above is the "Frequency Sampling Method" of filter design.

Comparison to "Windowing Method"

MattL in his answer has alternatively designed the filter using the "Windowing Method". (Both approaches use windowing, but the "Window Method" specifically windows samples of the true desired impulse response as the Inverse Discrete-Time Fourier Transform or IDTFT, while the "Frequency Sampling Method" windows the Inverse FFT derived impulse response). In general, when the "true" samples of the time domain impulse response are known, the Windowing Method may be able to compute the required filter with less resources and this is especially true if the target frequency response has abrupt steps. This is a trivial point when doing a one-time calculation of the filter coefficients using the resources of a personal computer, but can have significant consideration for applications designing dynamic filters on embedded platforms. To be clear, both will result in a filter with the same number of coefficients and similar performance with similar resource requirements to use the filter. This comment is just about resources to create the filter which in most applications is never a concern.

In summary for the two methods:

Frequency Sampling: Very simple. Sample the desired frequency response with a grid dense enough to minimize time domain aliasing, take the Inverse FFT to compute the estimated impulse response, truncate the resulting impulse response to the desired filter length and window. May require more computations to determine the resulting filter when many more Inverse FFT samples than final filter coefficients are needed. This approach may be preferred when we don't immediately know the time-domain impulse response and have no resource issues for designing the filter.

Windowing Method: (See MattL's answer) Opportunity for least computations to design the filter but requires knowledge of the actual impulse response. Simply window the impulse response if known at the desired length. This approach may be preferred when we know the time-domain impulse response, and when the resources or computation time to determine filter coefficients is actually a concern.

As an "apples to apples" comparison between the Frequency Sampling Method and Windowing Method for this specific case (the OP's spectrum), I recomputed the Frequency Sampling Method result using the same number of samples and the same Hanning window as MattL's Windowing Method answer (and important to reiterate, I did not for this comparison need to sample the frequency spectrum with a higher density of samples). The results suggest that, at least for this specific spectrum, the Frequency Sampling Method and Windowing Method are very similar in performance. This went against my past experience with getting significantly greater error without further oversampling the spectrum. I suspect this result and useful observation is specific to the smoothness of the OP's target spectrum and suggests the cases when Frequency Sampling can provide just as accurate of a result with no downside. ("Accurate" here means how well the resulting frequency response matches the desired frequency response).

The Python code used to create the filter coefficients using 'frequency sampling' is as follows:

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
import scipy.fft as fft
fs = 2      
N = 1001        # number of discrete freqs (using odd number for symmetry)
a = 1.2566      # constant
b = 0.00078956  # constant

k = np.arange(1, (N-1)/2+1)  # frequency indices for pos frequencies
f = fs/N * k                 # frequencies at sampled locations

#The spectrum
G = np.zeros(N)
hw = int((N-1)/2)   # half-way 

# create positive frequencies
G[1:hw+1] = a / np.sqrt(b*f[:hw])
# create negative frequencies as symmetric spectrum
G[hw+1:] = G[1:hw+1][::-1]
# DC Bin
G[0]= G[1]

# Create time domain response using IDFT
win = sig.kaiser(len(G),12) 
coeff = fft.ifftshift(fft.ifft(G)) * win    

# Recreate spectrum (as DTFT) and plot
w, h = sig.freqz(coeff, worN=2**16)
plt.figure()
plt.semilogx(f, 20*np.log10(G), label="True")
plt.semilogx(fs/2 * w/np.pi, 20*np.log10(np.abs(h)), label="Freq Sampling") 

Closing Thoughts

This post has presented the corrected approach to the OP's approach of using the Frequency Sampling Method for determining the time response from a frequency response, and provided a comparison to the alternative Windowing Method (detailed nicely in MattL's answer).

I have also gained through working through the OP's question a deeper understanding of the Frequency Sampling Method, in that there are cases (this one!) where the method provides just as good of a result as the alternative Windowing Method without any downsides. In most use cases for Frequency Sampling that I was previously familiar with, notably whenever the target spectrum has an abrupt discontinuity, the frequency spectrum must be oversampled to significantly more samples in order to achieve the same performance as the Windowing Method. As MattL has pointed out to me in the past, this does NOT mean the resulting filter need be longer, since the derived time domain response from the Frequency Sampling Method can be truncated and windowed. I was hoping to demonstrate that point in this post but found out surprisingly the two approaches have the same performance in this case without the need for increasing the density of frequency samples. My suspicion is the smoothness of the OP's target spectrum given the $1/\sqrt{f}$ shape.

Answer 2

As an alternative to the frequency sampling method discussed in Dan's answer, I will suggest another approach based on the analytical expression of the inverse discrete-time Fourier transform (IDTFT) of the desired frequency response. The resulting sequence can then be windowed and shifted such that the ideal response can be approximated by a linear-phase FIR filter. This is just the windowing method of filter design. Note that even though frequency sampling and windowing are usually not very popular for the design of frequency selective filters, they are very well suited for the approximation of smooth desired functions as the one given in the question.

In this answer it was mentioned that an analytical solution exists for the continuous Fourier transform, but that result cannot directly be used for a discrete-time implementation. Below I will derive the expression for the inverse discrete-time Fourier transform (IDTFT) of the desired response

$$G(\omega)=\frac{1}{\sqrt{|\omega|}}\tag{1}$$

Note that any desired scaling of $(1)$ can be directly applied to the IDTFT $g[n]$.

The IDTFT of $(1)$ is given by

$$g[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{e^{j\omega n}}{\sqrt{|\omega|}}d\omega=\frac{1}{\pi}\int_0^{\pi}\frac{\cos(\omega |n|)}{\sqrt{\omega}}d\omega$$

For $n\neq 0$ we obtain

$$\begin{align}g[n]&=\frac{1}{\pi \sqrt{|n|}}\int_0^{\pi |n|}\frac{\cos(x)}{\sqrt{x}}dx\\&=\frac{2}{\pi \sqrt{|n|}}\int_0^{\sqrt{\pi |n|}}\cos(x^2)dx\\&=\frac{2}{\sqrt{\pi}}\frac{C(\sqrt{\pi |n|})}{\sqrt{\pi |n|}},\qquad n\neq 0\tag{2}\end{align}$$

where $C(x)$ is the Fresnel integral $\int_0^x\cos(t^2)dt$.

For $n=0$ we obtain

$$g[0]=\frac{1}{\pi}\int_0^{\pi}\frac{d\omega}{\sqrt{\omega}}=\frac{2}{\sqrt{\pi}}\tag{3}$$

which also happens to be the limit of $(2)$ as $n\to 0$.

The IDTFT $g[n]$ is symmetrical and of infinite length. The corresponding filter is non-causal and unstable. In order to implement it approximately, we can apply a symmetrical window and then shift the response to the right such that the resulting filter is causal. Note that $g[n]$ decays only very slowly, so we have to use a rather long FIR filter to keep the approximation error small. However, if complexity and delay are of no concern, we can approximate the desired response as closely as we like.

In the example below I designed an FIR filter with $2001$ taps according to above procedure using a Hanning window. The figure shows the approximation error defined as the difference of the designed FIR response (in dB) minus the ideal response (in dB). Note that only the frequency region very close to DC is shown, because the error is virtually zero for frequencies further away from DC.

Below is the Matlab/Octave code to compute the filter coefficients of the FIR filter:

M = 1000;
N = 2*M + 1;    % number of taps = 2001
n = (1:M)';
% express (cosine) Fresnel integral by error function
x = sqrt( pi * n );
E1 = erf( ( 1 + 1i ) / sqrt(2) * x );
E2 = erf( ( 1 - 1i ) / sqrt(2) * x );
C = real( .25 * sqrt( pi / 2 ) * ( 1 - 1i ) * ( E1 + 1i * E2 ) );
g = 2 / sqrt(pi) * C ./ x;
g = [ g(M:-1:1); 2/sqrt(pi); g ];
w = hanning( N );
gw = w .* g;

Answer 3

If you want to get a real-valued $g(t)$, then its Fourier transform (FT) $G(ω)$ must be conjugate-symmetrical. So your expression of $G(ω)$ should be $G(ω)=\frac{a}{\sqrt{b|ω|}}$ after you set $\phi(ω)=0$. Otherwise, if $G(ω)=\frac{a}{\sqrt{bω}}$, it is impossible to get a real-valued $g(t)$ as its inverse fourier transform (IFT).

The IFT of $\frac{1}{\sqrt{|ω|}}$ has an analytic expression. That is,

$$IFT\{{\frac{1}{\sqrt{|ω|}}}\}=\frac{1}{\sqrt{|t|}}$$

You can use this result to directly get the expression of $g(t)$ in your case. Note that this is based on the definition of FT and IFT as

$$X(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}x(t)e^{-j\omega t} \, \mathrm{d}t,$$

$$x(t)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}X(\omega)e^{j\omega t} \, \mathrm{d}\omega.$$

If you are using a slightly different definition regarding the constants before the integral signs, there is a minor scaling adjustment.

Otherwise if you like, you may use ifft in MATLAB to obtain $g(t)$ numerically. You have to make sure to create a symmetric sequence of $G(\omega)$ before doing ifft in order to get correct result.