Is there a good method to extract a reference sin/cos from a group of signals? For the removal of a pure interference from a group of sampled signals, if I know the exact frequency of the interference I can simply use an I/Q demodulation technique (i.e., a phasor projection) to calculate it and remove it from the signals of interest. But I don’t have a pure version of the interference so I need to derive it from the data.

I know this is related to clock recovery, but with limited signal record lengths  any control loop (e.g., PLL or DLL) would be problematic.

So far I have simply hacked together an ad-hoc relatively long sequence of bandpass filtering, limiting, bandpass filtering, normalization for the I reference and afterwards all-pass filtering/Hilbert transform and normalization for the Q reference. It seems to work, but it’s too much of a hack. And it gets complicated if there are multiple interferers.

Is there a better/simpler way to extract these references?

If not, what would be a good alternative for simple bandpass filters? (as the requirements for these are quite stringent.)

If this was implemented in a real-time system (continuous data stream) instead, would you go the PLL/DLL route, or would this technique be good enough or even equivalent?

Although the most common method to remove these types of interference is to use a notch filter, the ringing in the impulse response of such filters is unacceptable for many applications.

This is particularly true for applications in which an intentional impulsive large artifact is present.

The best way to accomplish your desired result is to estimate the parameters of the interfering tone(s) and remove it(them) from the signal in the time domain.

To get you started on how to go about this, I recommend you read these two articles of mine:

Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT
Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

This needs to be done frame by frame, so there are a whole bunch of implementation details that come into play.  I'm really busy right now, so I will come back later with a fuller explanation.

For the two bin frequency formula, for best results, you will want to define your frame length on an integer + one half number of cycles.  That way you only need to calculate two DFT bins.  The amplitude/phase method will work on just two bins, but you will probably want to do four for that.

Given what you have said, this approach should be very plausible.

I would state that fact as you have 250 samples per cycle.  So suppose you aimed for 2 1/2 cycles, that would mean a frame length of 625 samples.  Any size in that range will put the interfering bin solidly between bin 2 and 3 (zero based indexing) and those are the only DFT bins you need so calculated a full DFT via FFT just for two bins is inefficient.  The longer duration you have, the more bin separation you get between frequencies, but less well the DFT works unless the tones are steady.  All the DFT cares about is sample counts, not what the underlying time scale is.

The proper units for frequency is cycles per frame which corresponds to the bin index or cycles per sample which you have given as 1/250.

Sample code is coming, and so is a further explanation.  You need to think in terms of frames, constructing a copy of your interfering tone in a separate set of buffer because you will need to have overlapped buffers.  Further details also forthcoming.

Still working on my two bin coding.  Hope you are too.

You want to think of your signal as a series of frames.  For each frame you are going to take a reading of the interfering tone.  Some frames will be good reads and some frames will be bad reads, and it is important to have a metric which tells you the quality of each read so you can set a threshold.

You will read the frequency and the amplitude and phase are best done as a cartesian complex number.  You will want to interpolate these values over a span of the good frames to give you parameters to generate a copy of the interfering signal so you can subtract it from your signal.

The reconstruction can be done by overlapping frames with averaging in the time domain as I said before, or a better way is to do a higher degree approximation of the angle in the sinusoidal function.

This method should work quite well for a slowly varying pure tone.  (and possibly the basis of a music compression algorithm?)


Pitch fluctuates unexpectedly in simple sine function I wrote a very basic C program to play a sine of a user-specified frequency. In the interest of portability, I have it spit values directly to stdout, so hopefully you can reproduce my problem on your own machine. This is the code I wrote:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

#define RATE 48000
#define TWOPI (3.14159 * 2.0)

int main(int argc, char** argv) {
        float frequency = atof(argv[1]);
        float phase = 0.0;
        float increment = TWOPI/RATE;
        int16_t sample;
        while (1) {
                sample = sin(phase) * 32767;
                putchar(sample & 255);
                putchar(sample >> 8);
                phase += increment * frequency;
        }
        return 0;
}

It can be compiled with:

cc sine.c -lm

and played with

./a.out 440.0 | sox -r 48000 -c 1 -t s16 - -d

That is: 16 bit mono audio at a sampling rate of 48,000 Hz.

I've tested this on both OpenBSD and MacOS and have noticed the same behavior: at around 6 seconds, the pitch of the wave jumps slightly. This happens again around 12 seconds. Further changes will take place as the wave continues to play. I have also replaced the sin() while loop with a prebuilt wavetable of a sine accessed through linear interpolation, only to notice the same problems.

Any ideas what might be off here?

You are using 32-bit floats, and not resetting the phase by subtracting 2 pi.  That means the phase will eventually climb up to where the sin() function's phase unwrapping algorithm doesn't have enough bits of quantization, or valid mantissa, left over afterward unwrapping.

Add this:

if (phase > M_PI) { phase -= 2.0 * M_PI; }

inside your loop somewhere.


Questions on the Generalized Tikhonov Regularization My first question is about the quadratic functional that is used in generalized Tikhonov regularization:
$$\Psi(f)=\frac{1}{2}\|f\|^2_\Gamma=f^T\Gamma f.$$
In the above equation what does $\Gamma$ represent, some special matrix?
Then, I don't understand why $\Psi(f)=\frac{1}{2}\|f'\|^2$ is preferable over $\Psi(f)=\frac{1}{2}\|f\|^2$ because
it imposes a penalty on the oscillations in the solution directly, rather than just a penalty on
the magnitude of the solution. Does that mean that the process is faster in the case we use derivative?
And one more question, I don't understand how to read the following figure:

Is the part on the right a result of denoising? Under the figure is written that we the covariance $C=\Gamma^{-1}=I$$ is used.

One way to interpret the Tikhonov Regularization is using the Maximum A Posteriori (MAP) framework.
Lets' say we have a model of the form:
$$ \boldsymbol{y} = H \boldsymbol{x} + \boldsymbol{n} $$
Where $ \boldsymbol{n} \sim N \left( 0, {\sigma}_{n}^{2} \right) $, namely Additive White Gaussian Noise, and the prior knowledge about $ \boldsymbol{x} $ is $ \boldsymbol{x} \sim N \left( 0, {\sigma}_{x}^{2} \right) $.
Then the MAP estimator is given by:
$$ \arg \min_{\boldsymbol{x}} \frac{1}{2 {\sigma}_{n}^{2}} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{1}{2 {\sigma}_{x}^{2}} {\left\| x \right\|}_{2}^{2} $$
So basically Tikhonov Regularization arises when we have prior knowledge or can model the parameters to have some kind of Gaussian Distribution.
Sometimes, our model is about the derivative (Very popular in Imaging). So we can say the derivative is distributed by Gaussian Distribution. Another way is having different covariance matrix (As in your example).


Is this signal perfect reconstructable? The question is as follow:

Let me do the analysis:

The downroad first:
\begin{align}
X_1(z) &= z^{-1}X(z)\\
X_2(z) &= \frac{1}{2}\left\{X_1\left(z^\frac{1}{2}\right)+X_1\left(-z^\frac{1}{2}\right)\right\}\\
X_3(z) &= X_2(z^2) = \frac{1}{2}\left\{X_1(z)+X_1(-z)\} = \frac{1}{2}\{z^{-1}X(z)- z^{-1}X(-z)\right\}\\
X_4(z) &= zX_3(z) = \frac{1}{2}\left\{X(z)- X(-z)\right\}
\end{align}
And the uproad:
    $$ X_6(z) = X_5(z^2) = \frac{1}{2}\{X(z) + X(-z) \}  $$

So the reconstructed signal is:
    $$ \bar{X}(z) = X_4(z) + X_6(z)  = X(z) $$

But I am very confused that this filter bank does not satisfy the perfect reconstruction according to PP V.'s Multirate systems and filter banks at page 193:

  $$F_0(z) =H_1(-z), \quad F_1(z) = -H_0(-z)\tag{5.1.7}$$

So I am very confused. Who could answer my question?

Yes the signal is perfectly reconstructed. Consider the process at each stage as I show using the block diagram below:

Consider each sample of the signal at each node in the diagram (each sample is shown using the sample index at the node for each row):

(Note: You see the same form of reconstruction in the FFT algorithm).

I will attempt to illustrate how the aliasing is cancelled (addressing MBaz's comment):

First consider the frequency response of $z^{-1}$ (and $z^{+1}$).  $z^{-1}$ is the z-transform of a DELAY of 1 sample in the time domain (see How/why are the $\mathcal Z$-transform and unit delays related?). The frequency response of a delay is constant magnitude, and linear phase (such as we would see with a cable of fixed length; very low frequencies would only be delayed a fraction of a cycle, while higher frequencies will be delayed several cycles:

Specifically for $z^{-1}$, since the delay is 1 sample per cycle, the amplitude will be constant (1) and the phase will be linear from 0 to $-2\pi$ for the frequency axis going from 0 to our sampling rate. (And $z^{+1}$ which is an advance instead of delay will be similar with positive phase versus frequency:

Phase only shown below (the magnitude is 1 for all frequencies)

Consider the following digital spectrum to see what happens as we pass through the system presented by the OP. (Note, when dealing with multi-rate systems, as well as mixed signal analog-digital systems, it has helped me to view the frequency axis from $-\infty$ to $+\infty$, which I depict below).

Now observe the beautiful thing that happens when we compare our original spectrum to that after a unit delay. I am emphasizing the that phase at $F_s/2$ has rotated 180°, and then fully 360° which is back to 0° at $F_s$.  This rotate is continuous from the linear phase shift from 0 to the sampling rate (which I have trouble showing in this picture-- I would be eternally grateful if someone could make me a 3d plot showing the actual spiral effect with the red and blue spectral artifacts since I use this plot in my class and as is takes too long to explain)

Now review what happens to the spectrum when we decimate by 2 (the aliasing that MBaz was astutely concerned about). For me, the analogy of the aliasing that occurs when you do analog to digital conversion has helped me intuitively see this process (since decimation is a "digital to digital conversion". That comment was just for those already familiar with aliasing in the A/D process.):

Finally, review what happens when we do a "Zero-Insert" interpolation. Zero-Insert interpolation perfectly maintains the spectrum with no distortion other than the image which now becomes part of our primary digital spectrum. This is because a zero-insert is done by convolving our waveform with the unit-sample response (impulse response).

With that important background, let's follow the system given by the OP. First we decimate our input signal into even and odd samples. Note the rotation caused by the delay and then the aliasing caused by the decimation (as described above, to really visualize this properly the spectrum after the delay would be a spiral where the blue is exactly 180° out of phase at the $f_s/2$ point only):

The even and odd samples are then interpolated by two by doing a zero-insert, and the bottom path is advanced one sample (so note that doing that is not a causal process, we would end up with a "parasitic delay" of one sample to be able to actually implement this, following what was done in the chart). Note that the $z^{+1}$ creates a spiral in the opposite direction. The combined effect of [$z^{-1}$, decimate by 2, interpolate by 2, $z^{+1}$] will shift the phase 180° for every frequency component above $fs/2$ in our original spectrum, while passing every frequency component below $fs/2$ with 0° phase shift (as well as create the images in the center of the new spectrum).  Thus by combining (adding) the two paths, the original spectrum is recovered!

Note: You can also subtract the two spectrums to get a reversed spectrum! Notice that the adding and subtracting is exactly what you do with a 2pt DFT: $$\begin{bmatrix}1 & 1\\1 & -1\end{bmatrix}\begin{bmatrix}x_1\\x_2\end{bmatrix}=\begin{bmatrix}y_1\\y_2\end{bmatrix}$$

We learn that aliasing can catastrophically ruin our spectrum and think that once the spectrums are contaminated we cannot separate the interference. That would certainly be the case here if we only had one of either path, but what is occurring is that our second path is keeping track of our aliasing in such a way that as long as we have both paths, we can constructively cancel out the interference.  This also gives insight into Multi-User Detection schemes.


Matlab freqz and custom implementation differences I am making some homemade tools on Matlab. I made function that plot the frequency response of a discrete transfer function (just like freqz does). When I compare my result with the result from freqz, I get something very similar, but not exactly the same.

Can somebody see why ? I am a little confused here.

Observations

I first thought it was related to frequency warping, but when I
applied warping on my frequencies, I got a similar result.
We also see that the phase does not wrap correctly, it passes from -179.3 to 180.7 (instead of a perfect -180 - 180);

Graph

Code

function [h,f] = py_freqz(b,a,fs,n)
%make sur both polynomials are same size
if (length(a) > length(b))
    b = [b, zeros(1,max(0,length(a)-1))];
elseif (length(b) > length(a))
    a = [a, zeros(1,max(0,length(b)-1))];
end
p = roots(a);
z = roots(b);

precision = fs/n;
w = linspace(0,pi-precision/2,n);
f = w/pi*fs/2;
h = ones(1,n);
phases = zeros(1,n);
unitcircle = exp(1i*w);
%Zeros
for k=1:length(z)
    v = unitcircle-z(k);
    h = h .* abs(v);
    phases = phases + angle(v);
end
%Poles
for k=1:length(p)
    v = unitcircle-p(k);
    h = h ./ abs(v);
    phases = phases - angle(v);
end

h = h * b(1);   % Input Gain

subplot(2,1,1);
xlabel('Frequency [Hz]');
ylabel('Magnitude');
plot(f,10*log(h));
grid on

subplot(2,1,2);
xlabel('Frequency [Hz]');
ylabel('Phase');
plot(f,phases/pi*180);
grid on
end

I generated the above graph using

py_freqz([0.09247185865855016 2*0.09247185865855016 0.09247185865855016 ], [1 -1.5668683171334845 0.9367557517676852], 1000,1000);

Several corrections:

This code does not make sense:

precision = fs/n;
w = linspace(0,pi-precision/2,n);
f = w/pi*fs/2;

Your precision should be the resolution in frequency domain (I think). But, now consider n=fs. Then precision=1 (which is fine), but w ranges only from (0...pi-0.5), i.e. you do not span the entire unit circle. The corrected code is (you actually dont need the precision variable):

w = linspace(0, pi, n+1); w = w(1:end-1);
f = w/pi*fs/2;

The trick is to create n+1 points in the interval [0,pi] and the remove the last one. This is much simpler than calculating the offset on your own.

You cannot add up the phases from each frequency zero and pole directly. Instead, you would need to do a modulo-addition (i.e. wrap around the phase). However, I suggest to simply track the phase+magnitude in a single variable h_all:

    unitcircle = exp(1i*w);
    h_all = b(1)*ones(1,n);
    %Zeros
    for k=1:length(z)
        v = unitcircle-z(k);
        h_all = h_all .* v;
    end
    %Poles
    for k=1:length(p)
        v = unitcircle-p(k);
        h_all = h_all ./ v;
    end

And finally you use this to extract phase and magnitude.

h = abs(h_all);
phases = angle(h_all);

Here's the full code for the function:

function [h_all,f] = py_freqz(b,a,fs,n)
%make sur both polynomials are same size
if (length(a) > length(b))
    b = [b, zeros(1,max(0,length(a)-1))];
elseif (length(b) > length(a))
    a = [a, zeros(1,max(0,length(b)-1))];
end
p = roots(a);
z = roots(b);

w = linspace(0, pi, n+1); w = w(1:end-1);
f = w/pi*fs/2;
unitcircle = exp(1i*w);
h_all = b(1)*ones(1,n);
%Zeros
for k=1:length(z)
    v = unitcircle-z(k);
    h_all = h_all .* v;
end
%Poles
for k=1:length(p)
    v = unitcircle-p(k);
    h_all = h_all ./ v;
end

subplot(2,1,1);
xlabel('Frequency [Hz]');
ylabel('Magnitude');
plot(f,10*log(abs(h_all)));
grid on

subplot(2,1,2);
xlabel('Frequency [Hz]');
ylabel('Phase');
plot(f,angle(h_all));
grid on
end

And the test script

b = [0.09247185865855016 2*0.09247185865855016 0.09247185865855016 ];
a = [1 -1.5668683171334845 0.9367557517676852];

[h, w] = freqz(b, a);
[h2_all, w2] = py_freqz(b, a, 1,8192);

%close all;
figure(1)
subplot(211);
plot(w/(2*pi), 20*log10(abs(h)), 'b-x');
hold on;
plot(w2, 20*log10(abs(h2_all)), 'r-');
hold off;

subplot(212);
hold off;
plot(w/(2*pi), angle(h), 'b-x');
hold on;
plot(w2, angle(h2_all), 'r-');

And the program output:


Regarding the choice of cost function in adaptive control - squared error vs absolute error I did search the question database regarding this question, and although one or two questions came close, they didn't really address my specific question.

In adaptive control based on minimizing tracking error (e.g. between plant and model), the designer is free to chose a cost function. More often than not the cost is selected as a function of the squared error.

But I've found in some practical applications that I can achieve a more robust controller by using absolute error. I understand that absolute error provides a more uniform weighting on the size of the error, and I suspect the squared error tends to 'wind-up' the adaptive controller with initially large errors. But I'm not sure how to show this in a generalized way. So I have two questions regarding this:

Is there perhaps a simple analysis that can demonstrate the stability characteristics between absolute and squared error choices in the cost function?
Any references on the matter?

This is an interesting question since both squared error and absolute error are convex functions, so they are both going to give the optimal solution when minimized. My intuition is that the $\ell_2$-norm (sum of squared values of error) converges to zero more quickly than the $\ell_1$-norm (sum of absolute values of error) when the search direction is right. For example, if the error is $\epsilon = y - Ax$. If $y$ is very close to $Ax$ and $0 < |y-Ax| \leq 1$, then $\epsilon^2$ is much smaller than $|\epsilon|$. Similarly, if $|y-Ax| > 1$, then $\epsilon^2 > |\epsilon|$. So the $\ell_2$ norm penalizes large errors more, and small errors less than the $\ell_1$ norm. This can also be understood by looking at the graphs of the two functions when $\epsilon \in \mathbb{R}$.

Of course the continuity of the derivative is another factor. The derivative of the absolute value technically does not exist at the minimum.


wavelet decomposition for time series signal Is it possible to use stationary wavelet decomposition as a tool to extract wavelet features for a time series? I can see how it works for image cases, but for a time series prediction problem say $x_i$, $i = 1,2,3,\cdots,n$ and you want to predict $x_{n+1}$, is it possible to use stationary wavelet decomposition to extract features and if yes why it makes sense?

A signal is just a 1D image. So if you can make it work for images, why not for signals? Pun aside, discrete wavelets combine multiscale smoother and differential operators, so they are been used as trend and singularity detectors for a while. Wavelet coefficients can be turned into feature vectors, directly or by condensating them across subband. For instance, wavelet subband coefficients are often believed to follow generalized Laplace-Gauss model distributions: $$A_{\alpha,\beta}\exp{\left(-\frac{|x-m|}{\beta}\right)^\alpha}\,.$$

The following coefficient histogram is for instance obtained from a stationary wavelet decomposition on a seismic signal (source: A Primal-Dual Proximal Algorithm for Sparse
Template-Based Adaptive Filtering: Application to
Seismic Multiple Removal)

And you can form a vector of subband features with estimated $\alpha_j $ and  $\beta_j $ at each scale $j$ for classification. Using stationary wavelets adds some shift invariance (see Difference between “equivariant to translation” and “invariant to translation”). Then, wavelet coefficients can be used as time-frequency anchors as in the Shazam music matching technique (An Industrial-Strength Audio Search Algorithm), or the deep-learning mimicking technique from the scattering wavelet transform developed by S. Mallat and coauthors.

Why does this make sense with (stationary)  wavelets?

Wavelets may act as multiscale change detectors, with a potential derivative behavior of wavelet functions at different scales, if you think about Haar (a discrete derivative) or the Laplacian of a Gaussian function with the Mexican hat wavelet,
They however can capture trends through the scaling function, since mother wavelets have vanishing moments,
They somehow can separate structured data from some noises, with orthogonality or whitening properties,
With stationary wavelets, the related features become integer shift invariant.


Is there a special reason why the principal range is defined as $-0.5\le F \le 0.5$? According to Digital Signal Processing by Ashok Ambardar, The principal range is defined as follows:

  Put another way, a DT sinusoid is periodic in frequency (has a periodic spectrum) with unit period. The range $-0.5\le F \le 0.5$ defines the principal period or principal range.

I see that the major rationale behind this definition is the unit period. If only for the purpose of uniqueness, it seems to me that any interval with a unit period surely meets the requirements. I also think that $0\le F \le1$ is more fit for the purpose, because it is more easy to remember and aesthetically beautiful. In conclusion, my question is:

Is there a special reason why the principal range is defined that way?

The $1$-periodicity indeed guarantees that any unit interval will contain all the information. Beauty can be subjective, but remember that art and science in the Western world have been strongly influenced by the Greeks. In Polyclitus's Canon and the Idea of Symmetria, you can see that beauty secrets:

  lay in the mastery of symmetria

With the $[-0.5,\,0.5]$ interval, you only have one number to remember, and the DC or $0$ frequency in the middle better underlines symmetry in quantities related to spectra.

Actually, beauty and simplicity put aside, the notion of "principal range" could be inherited from complex analysis, where the principal range of the phase is generally defined as $]-\pi,\,\pi]$. This is however a convention: any particular choice, defined in advance, can be called a principal range.

The $[0 ,\,1]$ frequency axis, often used by default by some FFT programs, is troublesome for many DSP learners, at the beginning. However, the  $[0 ,\,1]$  is still used quite often.

Note that the interval can be (should) open at one end: $]-0.5,\,0.5]$ or $[0 ,\,1[$.


Advantage of applying a Window Function in Analysis and Synthesis of STFT? I'm studying Short-time Fourier Transform and I learned that if you apply a windowing other than rectangular, you need to overlap the frames---for example, 50% overlap for Hamming family---to reconstruct the signal properly (like add them up to 1). However, I only need to apply this window function during the Analysis part and I can reconstruct the signal properly during synthesis.

But in some literature (like Audacity's Noise Reduction) I saw some methods that apply a second window function during the synthesis part and employ a 75% overlap-add. What's the advantage of doing the window function twice (analysis/synthesis) with 75% overlap compared to a single window function (analysis) with 50% overlap? I can't seem to find any explanation justifying the advantage between these two windowing techniques.

Edit: I also tried to implement the 75% overlap with Hann window applied on both analysis and synthesis part and I was able to yield a x1.5 of the original signal's amplitude. Did I do it correctly?

The STFT is an instance of an Analysis/Synthesis system. If we restrict to windows of finite length, it amounts to windowing blocks of length $L$, taking the DFT, and hopping by $K$ samples to perform the same operation. An A/S system is said invertible, or perfect, if you can recover the original signal from the blocks of DFT coefficients.
If $K>L$, some samples are missed, there are holes between blocks, and generally some information is lost (unless you have additional information).
If $K=L$, blocks are contiguous, and do not overlap.  Here, the transformation is critical, or non-redundant: you have the same numbers of Fourier coefficients than of samples. In that case, you can invert it even with non-rectangular windows. It suffices that the windows don't vanish, because windows without zero-values can be inverted (and Fourier as well).  There is only one inverse in that case. This is for your first question.
If you allow  overlap $K<L$, then several things may happen. I shall come back later, but basically:

a proper synthesis window can make a system with a vanishing analysis window invertible;
it can reduce artifacts resulting from  performing processing after the analysis stage (frequency removal, thresholding etc.);

In the continuous time domain, if you have an analyzing window $h(t)$, for any synthesis window such that
$$ \int_{-\infty}^{+\infty}g(t)h^*(t)\mathrm{d}t = 1$$
then the continuous STFT is invertible. Several observations:

These are very mild conditions to satisfy. This should not be surprizing, the STFT is extremely redundant, so the information is everywhere in the time-frequency plane
yes, you need a second window in the synthesis time-frequency plane (but it can be the same)
The windows $h$ and $g$ are complementary in second-order. Many possible $g$ exists, and especially:
$g=h$, provided that $ \int_{-\infty}^{+\infty}h(t)h^*(t)\mathrm{d}t = 1$ which is easy if the window does not have zero energy.

In the discrete-time/discrete-frequency domain, this is far more complicated in the general case. Yet, for power-of-two frame length, for certain window (triangle, Hann) and overlap ($1/4$, $1/2$, $3/4$), there exist direct inverses, easily computed from the analysis window. One shall compensate a bit in amplitude (like in your case), by a factor related to overlap and redundancy. Those cases are often standard, as they are quite easy.
However, there are generalized inverses, that allow for less redundant decompositions, at a cheap prices.
Related answers:

Meaning of Hop-Size in Filter Bank Interpretation of Short Term Fourier Transform
Normalizing STFT output and calculating energy of signal


Multidimensionnal wavelet admissibility criteria I've read there are two admissibility criteria for wavelets, both of which are designed to preserve total power of the signal (source: http://en.wikipedia.org/wiki/Wavelet#Mother_wavelet, as well as various scientific papers)

Condition for zero-mean: $$\int_{-\infty}^{\infty}\psi(t)dt=0 $$
Condition for square norm one: $$\int_{-\infty}^{\infty}|\psi(t)|^2dt=1,$$

where $\psi(t)$ is the wavelet kernel.
This leads to having a normalization factor of $\frac{1}{\sqrt{2\pi}}$ for the Gabor/Morlet wavelet.

Now, my question is: how do these admissibility criteria apply to 2D wavelets? (Let's say $\psi(x,y)$.)

My guess would be:

$$\int_{x=-\infty}^{\infty}\int_{y=-\infty}^{\infty}\psi(x,y)dx dy=0 $$
$$\int_{x=-\infty}^{\infty}\int_{y=-\infty}^{\infty}|\psi(x,y)|^2dxdy=1,$$

but I get inconsistent results in my computations.

Is my guess right? If not, what are they? Can you provide a reliable source?

It turns out the admissibility criterion does not dictate that the energy of the wavelet must be unity. To be classified as a wavelet, a wavelet must follow the following criteria:

The wavelet must have finite energy, so $E = \int_{-\infty}^{+\infty} |\psi(t)|^{2} dt < \infty $
The second condition is that the wavelet must have zero mean. (The intuition behind this is because the wavelet is acting like a matched filter, and we care only about how nicely the shape of the received signal matches the wavelet, and not the energy of the received signal). So if $\hat\psi(f) $ is the fourier transform of your wavelet, then: $$ C_g = \int_{0}^{\infty} \frac{|\hat\psi(f)|^2}{f} df < \infty $$ The $C_g$ here is called the admissibility constant, and the above property is called the admissibility criterion. This implies that $\hat\psi(0)=0$, because if it didn't, the above integral would blow up as you can see.
Finally, for complex wavelets, the last criterion states that the Fourier transform must be both real, and vanish for negative frequencies.

Now notice, no where does it state that the square norm of the wavelet must be equal to unity. It only states that the $E < \infty$.

In your case, it sounds as though the energy of your signal in the Wavelet domain, is not matching the energy of your signal in the Fourier domain. My best guess without further information is that you not normalize your wavelet to unit energy, since it is not demanded by the admissibility criterion. (First check though, if the energy of your signal in the spatial domain actually matches the energy of your signal in the wavelet domain). I would start there.


Overlap/Chunk processing possible with Zero-phase? I'm relatively new to DSP and I'm currently investigating FIR filter and IIR filters. From what I've found FIR filters can be implemented efficiently using the overlap-save method, but I was wondering if it might be possible to do the same using a IIR zero-phase filter?

For example say I have a fairly long signal if I divided the signal into overlapping chunks similarly to how overlap-save does it, and I zero-phase filter each of these chunks to prevent phase distortions then discard the overlapping samples as per the algorithm, is there any way I would be able to end up with a consistent end result? I've been looking around the web and in books but I haven't been able to find any sources that says it's possible except for this site:

http://www.dsprelated.com/freebooks/sasp/Overlap_Add_OLA_STFT_Processing.html

Where it says for Zero-phase to overlap the segment with two of it's neighbors. I've tried experimenting with the instructions and attempted to implemented and used an overlapping low-pass zero-phase filter on a signal as seen below, with the overlap length being 150 on each end of the segments. I've also padded the borders of each segment with 0s to prevent border distortions. From the looks of the result I got it doesn't seem like that it'll be smooth however when I combine them compared to if I filtered it all in one go.

If anyone could give me some pointers on what I might be doing wrong and where I need to look further it'd be great thanks!

If you apply a zero-phase IIR (for instance, by filtering both forwards and backwards) to a segment of a signal that has been zero-padded on both sides, and the two zero-padding lengths are equal or greater than the length of the IIR filter's impulse response (until that response drops below your desired noise floor), then you can sum non-overlapping segments, but you have to keep and include the zero-padded tails in the sum (similar to overlap add).  The error will be related to the (infinite length) ends of the IIR tails that have been cut off by finite length zero-padding.


How to equalize an ISI channel when the transmitted symbols are unknown? I am using MATLAB to simulate transmission of binary data sequence, $\mathbf{x}={x_k\in\{+1,-1\}\forall 0\leq k\leq N}$, over an ISI channel whose impulse response $\mathbf{h}=[h_0 h_1 \cdots h_{\mu}]$ is perfectly known. My channel observation sequence is $\mathbf{y}$ is a noisy version of filtered $\mathbf{x}$ that can be expressed mathematically as

$$y_k=Hx + w_k,$$

where $w_k\sim\mathcal{N}(\mu,\sigma_w^2)$.

For such system, equalizers can be designed to mitigate the effect of ISI. This often includes computing quantities like autocorrelation matrices $R_{xx}=E[x^{}x^T]$ and $R_{yy}=E[y^{}y^{T}]$, and cross correlation $R_{yx}=E[y^{}x^{T}]$

In MATLAB, I have access to $\mathbf{x}$ and I know noise variance $\sigma_w^2$ so computing $R_{xx},R_{yy},$ and $R_{yx}$  is straight forward.

My question is how do we compute such quantities in practice when we only have access to $\mathbf{y}$, our channel observation and  we don't know what input sequence $\mathbf{x}$ is?

The classic approach is to include a "training sequence" along with the data. The training sequence is known to both transmitter and receiver, so it can play the role of $\mathbf{x}$ in your question.

Now, in practice most wireless channels fluctuate over time. This means that the vector $\mathbf{h}$ is not fixed, but changes randomly in time. If the symbol rate is not too high, then the channel will remain essentially fixed during the transmission of several symbols. This time is called the "channel coherence time". When the symbol time is much smaller than the coherence time, we say the the channel exhibits "slow fading".

Over a slow fading channel, the training sequence is repeated approximately once every coherence time. During one coherence time interval, the equalizer remains fixed.


factoring poles / zeros: off by constant gain compared with textbook (From Schaum's DSP outline, 2nd edition, problem 5.32)

Book says factor it and extract H(z) from the factored product:

$$
H(z)H(z^{-1})= \frac{ \frac{5}{4} - \frac{1}{2}z - \frac{1}{2}z^{-1}  }{ \frac{10}{9}- \frac{1}{3}z - \frac{1}{3}z^-1 }
$$

ok.. no problem... so i convert to the form used by the book for nearly every problem:

$$
H(z)H(z^{-1})=  \left( \frac{z^{-1}}{z^{-1}} \right) \left( \frac{ \frac{5}{4} - \frac{1}{2}z - \frac{1}{2}z^{-1}  }{ \frac{10}{9}- \frac{1}{3}z - \frac{1}{3}z^-1 } \right)
$$

$$
H(z)H(z^{-1})= \left( \frac{ - \frac{1}{2}  + \frac{5}{4}z^{-1} - \frac{1}{2}z^{-2}  }{ \frac{1}{3}  + \frac{10}{9}z^{-1}-  \frac{1}{3}z^{-2} } \right)
$$

Then I apply quadractic formula to determine the zeros and poles:

$$zeros = \left\{ \frac{1}{2}, 2 \right\}$$
$$poles = \left\{ \frac{1}{3}, 3 \right\}$$

Then I rewrite $H(z)H^{-1}(z)$ in factored form:

$$
H(z)H(z^{-1})=\frac{\left(1-\frac{1}{2}z^{-1}\right)(1-2z^{-1})}{\left(1-\frac{1}{3}z^{-1}\right)(1-3z^{-1})}
$$

At this point I realize that I have a different answer from the book which says, the answer at this point should be:

$$
H(z)H(z^{-1})=\frac{\left(1-\frac{1}{2}z^{-1}\right)(1-\frac{1}{2}z)}{\left(1-\frac{1}{3}z^{-1}\right)\left(1-\frac{1}{3}z\right)}
$$

I think... well... ok.. it should be equivalent to my version... i can rewrite my equation to get it to match the book.. so i try that:

$$
H(z)H(z^{-1})=\frac{\left(1-\frac{1}{2}z^{-1}\right)(-2z^{-1})(1-\frac{1}{2}z)}{\left(1-\frac{1}{3}z^{-1}\right)(-3z^{-1})(1-\frac{1}{3}z)}
$$

$$
H(z)H(z^{-1})=\left(\frac{2}{3}\right)\frac{\left(1-\frac{1}{2}z^{-1}\right)(1-\frac{1}{2}z)}{\left(1-\frac{1}{3}z^{-1}\right)(1-\frac{1}{3}z)}
$$

At this point i'm scratching my head wondering why it doesn't match because its off by a different gain?  Wondering why the gain doesn't matter when picking out H(z) from the factorization of $H(z)H(z^{-1})$?  Is this good practice to write it the other way around? which way is correct?

Spectral factorization is not unique, unless you specify further constraints.

So let your given $H(z)H(z^-1)$ be (taking from your derivation) :

$$
\begin{align}
H(z)H(z^{-1})&= \frac{ \frac{5}{4} - \frac{1}{2}z - \frac{1}{2}z^{-1}  }{ \frac{10}{9}- \frac{1}{3}z - \frac{1}{3}z^-1 } \tag{1}\\
&= \frac{\left(1-\frac{1}{2}z^{-1}\right)(1-2z^{-1})}{\left(1-\frac{1}{3}z^{-1}\right)(1-3z^{-1})}
\end{align}
$$

Then all of the following are valid possibilities that satisf Eq.(1)

$$
H_1(z) =\frac{1-\frac{1}{2}z^{-1}}{1-\frac{1}{3}z^{-1}}
$$

$$
H_2(z) =\frac{1-2z^{-1}}{1-3z^{-1}}
$$

$$
H_3(z) =\frac{1-\frac{1}{2}z^{-1}}{1-3z^{-1}}
$$

$$
H_4(z) =\frac{1-2z^{-1}}{1-\frac{1}{3}z^{-1}}
$$

So you must have further constraints that make one of them as a unique answer. Such as minimum-phase, maximum phase, causal, stable, bothe etc...


Rough estimate of how close a human speaker is to a microphone? Is there a known technique for estimating how close a human speaker is to a single microphone? I'm hoping that there is a known technique that can analyze an incoming audio signal that contains a human voice with typical ambient room noise, and give a rough estimate of how close the speaker is to the microphone. Is there an algorithm/formula that can analyze such a signal and look for early reflections/reverb and then give a quality measurement back? Note, I am not looking for an estimate in units of distance (3 feet, etc.), just something that would tell me if the speaker is right on top of the mic, near the mic, far from the mic, or really far from the mic. Also, if it simplifies things, I can safely make the assumption that the speaker's voice will be significantly louder than any other human voice in the room.

Any papers or especially open source libraries would be fantastic.

PREREQUISITE: I would assume that an algorithm/formula that can detect if there is a human voice of significant volume present in an audio source would be a prerequisite for making the above quality measurement. Any information or resources along those lines would also be greatly appreciated.

Here is one possible suggestion: use a second order microphone (e.g. cardioid microphone). The human voice is a point source and has strong reactive velocity near-field component which velocity mics are sensitive too. In plain English this means you get substantially more relative bass when the speaker is close by. You could leverage to compare the energy in the 500-1kHz Band to the energy in the 100-200Hz band. With some reasonable calibration you ought be to able to tell whether the speaker is close or far.


Hexagonal Modulator Comparison with Square Modulator Hexagonal shaped Modulator is designed to get better BER than Square Modulator. The constellation diagram is like the figure:

It is supposed to give better BER measurement than Square as it has wider decision boundary. I have designed the modulator in this way.

Modulator for hexagonal

a2=[-2 -2 0 0 2 2 4 4];
a2=a2(:);
d=[1 -3 1 -3 1 -3 1 -3];
z1=a2+j*d(:);
z=[z;z1];
%plot(z,'O')

for i=1:length(s)
lodd=s(i);
ModDat(i)=z(lodd+1);

Demodulator for Hexagonal

a=[-3:2:3];
b=[3 -1 3 -1 3 -1 3 -1];
a=[a;-a];
a=a(:);
z=a+j*b(:);
%a2=[a2;-a2];

a2=[-2 -2 0 0 2 2 4 4];
a2=a2(:);
d=[1 -3 1 -3 1 -3 1 -3];
z1=a2+j*d(:);
z=[z;z1];
for fd=1:length(s)

for l=1:16; tes(l)=z(l)-s(fd); end
abs_dat_c=((abs(tes)).^2)';
min_dat=min(abs_dat_c);

for ds=1:16; if min_dat==abs_dat_c(ds); demoddat(fd)=ds-1;
    end;
end
demoddat=demoddat(:);

BER Measurement

N=170000;                       % No Of Symbols
M=16;
Sr=randint(N,1,[0,(M-1)]);      % random data form 0-15
s1=hexmod(Sr,M);  s1=s1(:);             % generate 16-qam data
% s2=qammod(Sr,M,0,'gray');s2=s2(:);
s2=qammod(Sr,M);s2=s2(:);

SNR=1:1:15
for k=1:length(SNR)
%% add noise---> LPN; CD; AWGN for ASE %%
Square_AWGN_Theory=awgn(s2,(SNR(k)+(10*log10(log2(M)))),'measured'); % Square Modulation
Hexagonal_AWGN_Theory= awgn(s1,(SNR(ppn)+(10*log10(log2(M)))),'measured'); % Hexagonal modulation

%% Demodulation Hexagonal
y1_Theory_hex=hexdemod(Hexagonal_AWGN_Theory,M);
y1_Theory_hex=y1_Theory_hex(:);
% [no_of_error_Theory_hex(ppn)

BER_Theory_hex(ppn)]=biterr(Sr(20000:length(y1_Theory_hex)-100),y1_Theory_hex(20000:
(length(y1_Theory_hex)-100)))
[no_of_error_Theory_hex(ppn)
BER_Theory_hex(ppn)]=symerr(Sr(20000:length(y1_Theory_hex)-100),
y1_Theory_hex(20000:(length(y1_Theory_hex)-100)))

%% Square Demodulation
y1_Theory_sq=qamdemod(Square_AWGN_Theory,M);
y1_Theory_sq=y1_Theory_sq(:);
% [no_of_error_Theory_sq(ppn)
BER_Theory_sq(ppn)]=biterr(Sr(20000:length(y1_Theory_sq)-100),
y1_Theory_sq(20000:(length(y1_Theory_sq)-100)))
[no_of_error_Theory_sq(ppn)

BER_Theory_sq(ppn)]=symerr(Sr(20000:length(y1_Theory_sq)-100),
y1_Theory_sq(20000:(length(y1_Theory_sq)-100)))
end

Output BER Comparison

Symbol error Rate Measurement

My Questions/Query

It is supposed that Hexagonal must give better result but i am not getting. although SER (symbol error rate is almost same ) but BER is not perfect. so please give me some suggestion by which i can get better result.

How to make differential coding for this Hexagonal Modulator ?

How to make the theoretical calculation to get a proper equation for this hexagonal modulator ?

For a rectangular QAM constellation, calculation of the symbol error rate is straightforward.
Exact calculation of the bit error rate is somewhat more complicated, but approximations
are more easily calculated based on the following assumptions.

Each signal point has four nearest neighbors (interior point)
or three nearest neighbors (edge points), or two nearest neighbors (corner points),
Note that 16-QAM has 4 interior points, 8 edge points, and 4 corner points.
The probability that a transmitted signal point is demodulated into
a nearest neighbor is orders of magnitude than the probability that it is
demodulated into a more distant point
A syymbol error results in exactly one bit error (that is, data bits
are assigned to signal points in Gray code order) so that the data bits assigned
to a signal point differ from the data bits assigned to a nearest neighbor
in just one bit position.

For the 16-ary hexagonal constellation with regular equiangular hexagons
(similar to but not quite the same as the constellation shown in the question),
the 4 interior points have 6 nearest
neighbors instead of 4 while the edge points have 4 nearest neighbors on average (some
have 3 and some have 5, others have 4 exactly), and two corner points have
3 nearest neighbors.  This has two effects. First, for a given spacing
between signal points, a symbol error is more likely since there are more
neighbors. Second, a symbol error can result in more than one bit error.
Thus, it is necessary to adjust the spacing between nearest neighbors
to compensate for this if one wants to have the same BER, and this will
typically increase the signal power needed.
In short, it is not immediately obvious that

  It is supposed to give better BER measurement than Square as it has
  wider decision boundary

is a believable claim.  In any case, the width of the decision
boundary is not quite as relevant to the BER calculation as the
distance of the constellation points from the decision boundaries,
and the number of nearest neighbors that each signal point has.

I don't decipher MATLAB very well, but if the picture with the
16 blobs was created by the MATLAB program included, then I
don't think that the constellation that has been created is
a (regular) hexagonal constellation at all. It looks for all
the world like you began with a rectangular 16-QAM signal constellation
(which has points at $(\pm m, \pm n)$ for $m, n \in \{1, 3\}$)
and moved the second and fourth rows by $1$ to the right, so
that the $x$ coordinate is $-2, 0, 2, 4$ instead of $-3, -1, 1, 3$.
Note that this increases the maximum signal amplitude in
one branch (I in this case) from $3$ to $4$ so that
the amplifiers have to have greater dynamic range and
the energy per bit also has increased. Any comparison
of the BER performance must take into account the fact
that the hexagonal constellation shown is operating at
a larger SNR than the standard rectangular 16-QAM constellation.
In a regular hexagonal constellation (which is what my remarks
above apply to), the points are at the
vertices (and centers) of regular (equiangular) hexagons and
are equidistant from their nearest neighbors, while
your distances are $2$ and $\sqrt{5}$ respectively.


Frequency response of each component of a system given its global response Given the following block diagram, find the frequency responses $H1(f)$ and $H2(f)$. The frequency response of the whole system has to be $H(f)=(\alpha_0+\alpha_1e^{-j2\pi T_1f}+\alpha_2e^{-j2\pi T_2f})^{-1}$

The fact that there's a loop confuses me. I would express $Y(f)$ as $Y(f)=X(f)H1(f)+(X(f)-X(f)H1(f)H2(f))\ H1(f)+...$, but it doesn't seem to be correct. Could you give me some hints? Thanks in advance!

  Solution: $H1(f)=\alpha_0^{-1},\ H2(f)=\alpha_0^{-1}\alpha_1e^{-j2\pi T_1f}+\alpha_0^{-1}\alpha_2e^{-j2\pi T_2f}$.

With feedback systems such as the one given it's often easy to define an additional signal at the output of the adder. This gives the following equations:

$$U(f)= X(f)-H_2(f)Y(f)\tag{1}$$

and

$$Y(f)=U(f)H_1(f)\tag{2}$$

Now you can solve Eqs $(1)$ and $(2)$ to get the frequency response $H(f)=Y(f)/X(f)$.


Why does part of my image turn black when applying functions to a part of an image in Matlab? When I select a part of an image and apply a negative function to it, it works 100%:

When trying to apply the same method using another function, the part of image becomes black:

What is the problem and how can I fix it?

'edited image' is of type uint8 and the selected area which you have taken for processing finally after applying all operations is coming out as type double.But when you place it in original image it is considered as of type uint8 because the image is of type uint8. Hence before putting it in original image you should multiply x by 255 so that it get be properly scaled.


How does an IIR system affect magnitude and phase of a sinusoidal signal Consider an IIR system with impulse response $h[n]=\left( \frac{1}{\sqrt{3}} \right)^n u[n]$. If I apply $x[n]=\cos(n \frac{\pi}{2} + \varphi)$ at the input, how can I determine the change in magnitude and phase of the output signal, i.e. what is $A$ and $\phi $ in $y[n]=A\cos(n \frac{\pi}{2} + \varphi + \phi)$?

HINT:

Compute the system's frequency response $$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}$$
Compute its magnitude and phase
Figure out (by searching this site or any textbook on Signals and Systems) how the magnitude and phase of $H(e^{j\omega})$ affect a sinusoidal input signal.


Decimating physics data At work we run a physics simulator (video game) that operates in an inner loop at 360 Hz.  Data from that inner loop is collected 60 times a second and sent out to hardware (a force feedback steering wheel) in order for playback to the user.

I have two questions on this system.

First what is the best way to decimate the 360 Hz data down to 60 Hz.  Right now we do a straight average of 6 samples together, but I suspect that this is causing some aliasing of the data and we are very sensitive to high frequency noise.  We are also extremely sensitive to latency and even adding in 4 samples of latency (in the 60 Hz loop) would cause a noticeable delay between the onscreen display and the force feedback device that the users would pick up on.

Second, the hardware (steering wheels) is typically made to very poor tolerances and usually has a difficult time playing back data above about 20 Hz before the signal significantly degrades.  Would it be better to low pass filter the 360 Hz data down to 20 Hz, or to low pass the 60 Hz data to 20 Hz, and what would be the most efficient way to handle this.

I currently have a single pole IIR filter that you can optionally define on the 60 Hz signal, this does smooth the data out, and adds in a considerable amount of latency as well.  Mixing less than 0.20 of the current sample causes an odd disconnected feeling in the wheel so that is my lower bound.  I would assume that going to a more complex filter would give me a much faster falloff and would both smooth out the high frequency chatter and noise in the wheel without increasing the latency beyond 4 samples (66 ms)

I'm half way through Steven W. Smiths book "The Scientists and Engineer's Guide to Digital Signal Processing" so I have a limited understanding of DSP's.  I do have 20+ years of programming experience in 2D image manipulation so I have been using dsp's all along without knowing what they were.  The point is go easy on my poor mind :)

I should also point out that the 60 Hz loop is running in 'real time' that is timed to the pc clock, and all communication with the wheel is also in real time, however the inner loop is processed 60 times a second and there is no way to get access to the inner loop data at a faster rate than 60 Hz.

This is a really very nice applied signal processing question. Unfortunately your lack of DSP and physics expertise, builds some barriers on the way of your understanding. Hopefully this answer helps to tunnel through.

Now, your first question is about downsampling a discrete time signal. In fact there are already plenty of such questions here under the title of interpolators, downsampling etc. These are already answered with great detail. You better make a search on the site for more detailed answers but here once more I'll give a brief reply.

In downsampling a discrete-time signal (say from 360 hz to 60 hz in your case, by a facttor 6) the most critical thing to consider is the Spectral Utilisation of the original (360 Hz) signal. Based on Shannon-Nyquist sampling theorem criteria, this signal is assumed to contain meaningful (non-trivial) information in the band of frequencies 0 to 180 hz. However for many applications the frequency band of interest within the say 0-180 hz range will be a much smaller band, say for example from 0-20 hz. In such a case, the original signal is said to be "oversampled" (by a factor of 9 in this case). When a signal is oversampled by a factor of say M, then it is somewhat acceptable to discard every (M-1) samples and use the Mth one to represent the information content of the original signal without loss (except at the higher end of the band which gets to limit of aliasing). But this will/may produce some noise due to aliasing being on the very edge of it.

There won't be any aliasing, as long as the out-of-interest frequency band in the original signal contains no energy. In order to guarantee this condition of zero energy in the out of interest band, it is required to filter out those out of interest bands by a suitable filter of sufficient quality, before decimating the signal.

The quality of such a filter depends on, therefore, i-the spectral distribution of the original signal, ii-band of interest, iii-downsampling factor and other requirements such as real time processing or offline, accuracy, CPU power etc.

In your case, the answer lies in the following: the best downsampler should produce the most accurate representation of the band of interest signals and reject any irrelevant information for your steering hardware to process.

The short answer is this: as long as your system is configured to use 60 hz information (30 hz band) for the force-feedback processing, it seems that (and from my experience with mechanical acceleration signals, it is indeed so) this band of signals is adequate to represent average mechanical control. (as long as I can remember, most mechanical acceleration sensors do use at most a band of a few hundred hz). So just use a simple low-pass filter of cutoff frequency at pi/6 for better performance.

The second questions simply asks which is better to choose i-filtering 360hz into 20hz or ii-60 hz into 20 hz. The answer is the latter. It will be a much shorter filter to apply, and the transition band would also be easier.

When the latency of a few samples becomes so critical as in your example of real time forced-feedback control of a physical simulator, then high accuracy filters becomes impractical to use. May be you should completely discard using any filters at all, as you have pointed out the physical limitations of the wheel mechanics already prevents bulding of responses faster than 20 hz...


How can I classifying this License Plate images? Matlab , LPR , I want to perform automatic character recognition on licence plates. One of my preprocessing step is binarization. However I have four different class of images and I have to use a different automatic thresolding algorithm on each class to get good binarization results.

I would like to choose the right algorithm automaticaly, my idea is to first perform classification on the licence plates images. In the following picture, each row represent a class.

picture:

http://www.mathworks.com/matlabcentral/answers/uploaded_files/31813/class.jpg

I'm working with MATLAB, how can I do that?? what is the best way?
Any ideas are appreciated.

Here's a procedure you should try in a first attempt :

You have RGB images, so each images is a 3*RC matrix, where R is the number of rows and C the number of columns.
We will consider each pixel of an image as a realization of a random vector $\mathbf{x}$ of dimention 3. Hence, having an image, you have RC realization of $\mathbf{x}$. With the classification exemple you gave, each classe has 4 images, so we have 4*R*C realization of $\mathbf{x}$. We are abble to estimate some stastical moments, that we hope will be good descriptor, such as mean, covariance etc ... (mean can be estimated with empirical mean or median).
Let's just stick with the mean for now. Let us call $m_1,m_2,m_3,m_4$ the means of each class you gave.
Now we have a new image we want to classify, we estimate it's 3-dimensional mean $\mathbf{\mu}$, then using nearest neighbour, the class is given by :
$$
C=\arg_i\min\left(\lVert\mathbf{\mu}-m_i \rVert _2\right)
$$

This is very basic.
You may improve this scheme using additional features like variance (diagonal of the covariance matrix) or covariance matrix. You may estimate variance or covariance using empirical variance or empirical covariance estimator (realy just the classical formulae). However the euclidian distance won't be a good distance since you will comparing inhomogenous vector (the features will be composed of mean + covariance), you'll want to use a more advanced distance such as Bhattacharyya distance for the special case of gaussian distribution.


PCMA signal separation I want to work on co-frequency and co-time blind signals separation or carrier-in-carrier signal separation project. In articles published under the heading of blind separation of these signals, it is always assumed that one of the signals is known and the other signal is obtained by using that. I couldn't find an article in which both signals were unknown and separated. This has confused me. Is there anyone who has worked on this and guided me? Is there anyone who has worked on this and guided me?

Patent articles indicate that carrier-in-carrier separation is done by using adaptive filter. In this way, the received signal (which is the sum of the first and second signals) enters the filter with one of them and by subtracting it from the received signal, another signal is obtained. Some articles have also pointed out that if the power difference between the two signals is significant, the second signal is obtained by using a trick.

However, various articles have suggested that if the power difference between the two signals is 0 dB, it is possible to obtain both signals blindly and without any of them. Is such a claim true?

This reference could be useful for blind equalization of CiC Signals


LMS and delta impulse response- Equalization concepts I have been reading about equalization and channel estimation in order to apply to my application of estimating trajectory path for a two-wheeled robot. The application is quite similar to equalization.  I am studying the Least Mean Square algorithm (LMS) and have some basic doubts regarding how to determine if equalization is proper or not i.e., if there is a measure to quantify how much of the signal has been estimated properly. I am following this MATLAB implementation.

Doubt 1 :
The cost function for LMS is constructed using the input, so we are minimizing the error function. So, based on my understanding the input is known to the receiver. Hence, the LMS equalizer is not blind. Please correct me if wrong.

Doubt2:
On running the code (using the same parameter settings in the code), I cannot understand what is the significance of the fourth tap having a magnitude one. Is this the delta impulse in the third subplot.

Does having a delta function mean that equalization is correct? I had read that for equalization to happen, the convolution of the estimated channel coefficients (here these are the weights of the LMS) and the known channel impulse should be a delta function. Does this graph convey the same? Could somebody please explain this concept of delta impulse and the graph?

You're right. LMS equalizer uses a known input to minimize the error. For communication purposes, this is either provided by a training sequence, or in a decision directed mode, the detector decisions are fed back as known data.

The delta function is also correct. Suppose that the channel impulse response is $h(t)$ and frequency response $H(f)$. Then the equalizer should ideally be the inverse of this function in frequency domain (remember the convolution in time domain and hence multiplication in frequency domain).
$$H(f)\cdot \frac{1}{H(f)} = 1$$
and its inverse Fourier transform is
$$F^{-1}(1) = \delta (t)$$

The fourth tap being 1 implies that there is a delay of 4 symbols with respect to your time base at the Rx. Finally, to further verify your results, plot the frequency responses of the channel and the equalizer response, and you will see an approximately inverse relationship.


Intuitively explain Bi-linear time frequency distributions, someone please? I have a fairly good understanding of Wavelet Analysis, but what are these bilinear distributions and how do they differ from Wavelet Transform?

As you mentioned - all these methods share common principle, they allow for representation of our signal in time-frequency domain. First thing to notice is that wavelets are very different from the STFT-like (linear) methods as they provide scale dependent frequency resolution (which in case of STFT is constant and governed by Heisenberg uncertainty principle). Probably you've seen this famous picture:

So called bi-linear transform mostly refers to Wigner Transform which can be considered as the Fourier Transform of signals' instantaneous autocorrelation function. It is governed by following equation:

$$W(t,f)=\int_{-\infty}^{\infty} x \left( t+ \frac{\tau}{2} \right) x^* \left( t-\frac{\tau}{2} \right) e^{-2\pi i f \tau} d\tau$$

where $^*$ denotes complex conjugate and $\tau$ is the time delay.

One might notice that in analogy to STFT, the shifted signal itself is considered as a window. It is also worth of mentioning, that Wigner Transform provides perfect time and frequency resolution at the same time - this is simply great for signals with fast changes of instantaneous frequency. On the other hand WT is the bi-linear transform, which implies major disadvantage: for each linear combination of two signals $x(t)$ and $y(t)$, result consists of two parts: auto-terms and cross-terms. This issue produces interferences in time-frequency representation known as "beats", which can make interpretation of WT difficult. If you are looking for more specific mathematical description, then please refer to article posted below - I see no point in rewriting equations.

Example is shown on a figure below, where two linear chirps were analysed. Cross-term is clearly visible in between, which makes interpretation of results troublesome. Sometimes it happens, that magnitude of cross-term is even higher than auto-terms.

Although some methods were described in literature how to deal with this problem. Please refer to the following article for more details.

  R.B. Pachori and P. Sircar - A new technique to reduce cross terms in the Wigner distribution


MIMO System Model in Matlab I have this lame doubt which has ended up confusing me which is probably trivial. I am modelling a MU-MIMO system with a Rayleigh fading channel. I have created a wanted signal as an OFDM signal, and generated a 3D Rayleigh Channel matrix over all sub carriers. I want to now multiply the transmit vector(OFDM signal) to the 3D channel matrix. Should i convert both the channel matrix and wanted signal in frequency domain and add noise

$y(f) = H(f)x(f)+n(f)$, and then convert it back to time domain
$\textrm{ifft}(y(f))$ and provide an estimate using MRC, ZF or MMSE methods.

$\hat{y}(t) = H(t)’ y(t)$

Or is it just a direct multiplication of

$y(t)=H(t)x(t)+n(t)$

as i understand from most of the literature (y = H x+n)is the usual system model in time domain. If so why isn't it convolved when its combined in time domain, similar to any communication system, where the input vector is convolved with the channel matrix to obtain the output vector in time.

In a wireless system, most of the noise is added by the receiver itself. The received signal is very weak (-100 dB or less is typical), and its power is comparable to the thermal noise of the receiver's amplifiers. That is the reason why noise is modeled as being added to the signal present in the antenna.

The noise still affects the estimated signal, of course; if the received vector is $r$, a naive ZF receiver would calculate

\begin{align} \hat{y} &= H^{-1}r \\ &= H^{-1}[Hy+n] \\ &= y + Hn. \end{align}

Since the noise $Hn$ is correlated, this receiver has terrible performance.


What is (Fourier frame length/2 + 1)? If $N$ is discrete Fourier transform's frame length, and $N/2$ is half the frame length, how would you call $\frac N2 + 1$?

frame_length      = N      ;
frame_length_half = N/2    ;
?                 = N/2 + 1;

Using your notations: number_frequency_bin, i.e. the number of "non-redundant" frequency bins, the first (corresponding to the DC) and last (corresponding to Nyquist) ones being always real for a real signal. The first is a (normalized) sum of samples, the last a sum of samples with alternated signs (since  $e^{i\pi n }=(-1)^n$).


How we calculate Precision-Recall Curve? As far as I know, precision and recall are two single values. How we can plot a curve from these two single values? I think I should calculate a set of values for each of them but how?

afterwards, the curve can be depicted by using the fact that they are inversely related.

$$
\text{Precision}=\frac{tp}{tp+fp}
$$
$$
\text{Recall}=\frac{tp}{tp+fn}
$$
Where $tp = \text{True Positives}$, $fp = \text{False Positives}$ and $fn = \text{False Negatives}$.

Anybody can explain how we can plot a curve from two single values?

One answer can be found here but I do not believe that I have caught the point well.

Obviously, you can't plot a curve from a single pair of values, so you have to find a parameter to vary that changes your problem in a way that makes sense to what you're trying to do.

For example, suppose your binary classification problem has a choosable threshold.  If you vary the threshold up and down, then it will change the precision and recall.

Then you might get a plot like the one below.  R code below that generated it.

R Code

#26568

Ndata <- 1000

c1 <- runif(Ndata,1.5,3)
c2 <- runif(Ndata,1,2.5)

data <- c(c1, c2)

thresholds <- seq(1,2.99,0.01)

precision <- thresholds*0
recall <- thresholds*0

for (k in 1:length(thresholds))
{
  threshold <- thresholds[k]
  cl <- data*0
  tp <- 0
  tn <- 0
  fp <- 0
  fn <- 0

  for (i in 1:length(data))
  {
    if ( (data[i] > threshold) && (i <= Ndata))
    {
      tp <- tp + 1
    }
    if ( (data[i] > threshold) && (i > Ndata))
    {
      fp <- fp + 1
    }
    if ( (data[i] <= threshold) && (i <= Ndata))
    {
      fn <- fn + 1
    }
    if ( (data[i] <= threshold) && (i > Ndata))
    {
      tn <- tn + 1
    }
  }

  precision[k] <- tp / (tp + fp)
  recall[k] <- tp / (tp  + fn)
}

plot(precision,recall)


What is spatial correlation and spatial convolution? What does the term spatial correlation and spatial convolution mean with respect to digital image processing ?

These terms exist mainly for historical reasons. In signal processing the signal is a one-dimensional function of time.  So people talk about the time domain vs. the frequency domain.  On the other hand, in image processing you are looking at a 2D function of $x$ and $y$, and there is no notion of time.  Instead your are talking about spatial frequencies.  Hence, spatial correlation and spatial convolution.

Typically, in image processing you simply talk about convolution and correlation.  The term spatial usually shows up when 2-D convolution and correlation are introduced to people with background in signal processing.


Virtual Microphone using Mic Impulse Response I am working on a virtual microphone project which tries to emulate the output of an audio signal that would be obtained using another microphone. Let me know if this sounds reasonable as a solution or not.

Let's suppose, we have two microphones: Mic1 and Mic2. A sound source s is played through a loudspeaker using Mic1. Let h1 and h2 be the impulse response of Mic1 and Mic2 obtained using the same loudspeaker respectively. Now, if we record using Mic1, we obtain:

s1(t) = h1(t) * s(t) % (*) operator is convolution.

If we want to achieve s2 (the sound produced by Mic2), we can do the following:

% Take fft of h1(t) & h2(t)
H1f = fft(h1); H2f = fft(h2);

% Obtain S(f) using division in fft
S1f = fft(s1);
Sf = S1f/H1f;

% Obtain s2(t) using multiplication in freq domain & ifft
S2f = Sf .* H2f;
s2 = ifft(S2f);

Does this sound feasible as a solution to the problem? A more practical implementation (suggestion) would be appreciated.

Will this work? Well, it will do something, but it will not magically make one microphone sound like another. A microphone is much more than a signal receiver with an impulse response. The character of a microphone is determined by its directional characteristic and how it couples to diffuse and direct field components, the kind of coupling to the pressure field, the nonlinearity of the membrane and integrated preamplifier, the noise of the membrane and the electronic components.

You will only be able to capture a small part of all that with an impulse response, and very likely exactly not the part of a microphone that is responsible for the praise of its sound.

But even if you ignore all that, you're still left with the problem of deconvolution that is not trivial, especially in the presence of noise of a potentially unknown character. Your spectral division approach is rather naive and will work to a certain degree. If you want to get good results, you will have to take a harder route and you will probably not even in that aspect get what you imagine. Also, you will always have to start with a very good microphone. A bad microphone does not allow for such a simple improvement usually.

In short, even if I sound overly pessimistic: No, this will not work how you probably expect it to. If you want to do this to start a research project then you're good. But it's not going to be a plug and play product for microphone aficionados.


First-order low pass filter with bounded derivative Given a function $x(t)$ and a first order low-pass filter $\frac{dy}{dt} = \frac{1}{\alpha} (x(t) - y(t))$, is there a way to choose $\alpha$ analytically such that the derivative of $y(t)$ is bounded, i.e. $\forall t \frac{dy}{dt} \leq c$ for some given constant $c$?

Thanks for your help.

Note that this question is similar to the one here.

I will show you two bounds on the derivative of the output signal, both based on the system's impulse response. First, by taking the Laplace transform of the time-domain equation and rearranging you get the Laplace transform of the output signal in terms of the transform of the input signal and of the system's transfer function:

$$Y(s)=X(s)H(s)=X(s)\frac{1/\alpha}{s+1/\alpha}$$

The inverse Laplace transform of $H(s)$ is the system's impulse response

$$h(t)=\frac{1}{\alpha}e^{-t/\alpha}u(t)$$

where $u(t)$ is the step function. In the following I will assume that $\alpha>0$ is satisfied, because otherwise the system is not stable.

The transform of the derivative $y^{\prime}(t)$ is

$$sY(s)=sX(s)H(s)=[sX(s)]H(s)=X(s)[sH(s)]\tag{1}$$

from which it can be seen that the derivative $y^{\prime}(t)$ can be obtained either from filtering the derivative of the input signal $x^{\prime}(t)$ with $h(t)$, or by filtering $x(t)$ with a system with transfer function $sH(s)$. The impulse response of this system is

$$h^{\prime}(t)=\frac{1}{\alpha}\delta(t)-\frac{1}{\alpha^2}e^{-t/\alpha}u(t)\tag{2}$$

From these two interpretations we can now derive two bounds on $|y^{\prime}(t)|$. First note that the magnitude of the output of a linear time-invariant system can be bounded as follows:

$$|y(t)|=\left|\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau\right|\le
\int_{-\infty}^{\infty}|h(\tau)||x(t-\tau)|d\tau\le
|x(t)|_{max}\int_{-\infty}^{\infty}|h(\tau)|d\tau\tag{3}$$

Since $y^{\prime}(t)$ can be obtained from filtering $x^{\prime}(t)$ with $h(t)$ we get the following bound

$$|y^{\prime}(t)|\le|x^{\prime}(t)|_{max}\int_{-\infty}^{\infty}|h(\tau)|d\tau=
|x^{\prime}(t)|_{max}\frac{1}{\alpha}\int_{0}^{\infty}e^{-\tau/\alpha}d\tau=|x^{\prime}(t)|_{max}\tag{4}$$

This is the first bound on $|y^{\prime}(t)|$, which is of course only useful if one knows an upper bound on $|x^{\prime}(t)|$.

The other bound is obtained by noting that $y^{\prime}(t)$ can also be interpreted as filtering $x(t)$ with the impulse response $h^{\prime}(t)$ given by (2):

$$y^{\prime}(t)=(x*h^{\prime})(t)=\frac{1}{\alpha}x(t)-\frac{1}{\alpha^2}\int_0^{\infty}e^{-\tau/\alpha}x(t-\tau)d\tau$$

This leads to the following bound

$$|y^{\prime}(t)|\le\frac{1}{\alpha}|x(t)|_{max}+\frac{1}{\alpha^2}|x(t)|_{max}\int_0^{\infty}e^{-\tau/\alpha}d\tau=\frac{2}{\alpha}|x(t)|_{max}\tag{5}$$

Combining (4) and (5) we obtain

$$|y^{\prime}(t)|\le\min\left\{|x^{\prime}(t)|_{max},\frac{2}{\alpha}|x(t)|_{max}\right\}$$

Obviously, both bounds depend on properties of the input signal $x(t)$, and only one of them depends on the system parameter $\alpha$. In sum, without knowing anything about the input signal $x(t)$, there is no way to choose $\alpha$ such that $y^{\prime}(t)$ is bounded for all input signals.


How to make sure that a sinusoid wave with different subsequent frequencies matches up? A very basic formula to transform a sequence of notes into a signal would be:

$$
x(t) = A\cos( 2  \pi  (2^{current\_note}/12) t)
$$
where $current\_note$ is $notes[\lfloor t / period\_per\_note\rfloor]$ and $a[k]$ is the indexing operation into a string; in this case into a sequence of numerical values representing notes on a chromatic scale.

Listen to an example

However, when making a signal like this, whenever the frequency changes (when moving from one note to the next), the function is discontinuous, often starting at a completely different amplitude because in general, the periods of the two sinusoids with different frequencies end up being very different (not being a multiple of the same 'base frequency')-

This can also clearly be heard in the example, as 'clicks'.

Now the question: How can we make sure that the sinusoid with the 'new' frequency starts at the same value of x as the one with the old frequency?

I am guessing that this somehow involves keeping track of the current 'phase', so that we can time-shift the new frequency to make sure their amplitudes match up. But I have no idea how to approach this, neither from a mathematical nor from a programmatic perspective.

You need to make sure that the phase of continuous.

The easiest way would be something like this.

x(n) = cos(phase);
phase = phase + 2*pi*current_frequency*sample_time;
current_frequency= update_frequency(n);
n = n + 1;

Note, that this is very inefficient code and only for illustrative purposes.


Ultrasonic TDoA / ToA Measurement? I am trying to make an ultrasonic flowmeter. The device has two ultrasonic tranducers (resonance frequency of about 200kHz) that sample one ultrasonic signal in the upstream directon and then immediately after one signal in the downstream direction. I know a correlation between the upstream and downstream signals might be enough but the problem is that we don't have a lot of processing power available. We have signals that are 1024 samples long each and a full cross correlation would take too long to process. The maximum sampling Rate is about 6MHz if that helps.
Since I don't come from an audio background, I am also finding it hard to figure out at what point the signal really 'arrives' and how to detect it at the resolution that I need. With the SNR that I have, it seems like the exact startng point would be lost in the noise?
No matter what I try, I seem to end up with too many wrong measurements where the time measured is exactly one wave period away from the correct value.
I know this question is far too vague but I would really appreciate any ideas on how I could understand the problem better!
Here's what the signals look like:

So the main questions are:

What is the exact time of arrival of these two signals at the ultrasonic transducer?
What is the differential time of arrival between both the signals?

The FFT-based cross-correlation is done as follows:
TX = np.fft.fft(tx)
RX = np.fft.fft(rx)
CORR = np.multiply(TX, np.conjugate(RX))
corr = np.real(np.fft.ifft(corr))

Here is a NumPy script that tries a solution.
It does the following:

create a TX and delayed RX signal that similar to your waveform (a sine wave modulated by an envelope cosine with zeros at the front and the back and noise superimposed
calculate the cross-correlation with SciPy's signal.correlate
calculate the cross-correlation with FFTs
find and print out the location of maximum correlation
if all went well, the resulting number should be the same as the delay between TX and RX

In the graph below, you see the time domain test signal at the top, and the correlation result at the bottom (blue = signal.correlate, orange = FFT-based).

As you can see, there's a shift between the blue and the orange result. Figuring out why is an exercise for the reader (I'm learning this stuff myself). It has to do with cyclical convolutions, zero padding etc. You should be able to get a better understanding than what I'm able to offer by googling "FFT convolution". Here is one such Google result that talks more about this.
In my script, instead of doing it fundamentally right, I went the pragmatic road and changed some fixed constant to make the convolution and the FFT convolution come up with the same result.
The full script finds the maximum value of the 2 correlations, like this:
    max_corr = np.argmax(corr)-len(tx)//2
    max_corr_fft = np.argmax(corr_fft) - len(corr_fft)

And prints the result:
max corr: -129
max corr fft: -129

129 matches the delay that was used to generate the RX signal.


Kalman Filter for position : introducing acceleration estimates I want to estimate the position on a 3D environment by introducing only acceleration estimates. Is that possible? If I use the extended Kalman Filter and introduce these estimates will I have the appropriate results?
Thank you in advance!

If the acceleration readings are performed in an inertial (non rotating in particular) frame then yes it's sufficient to use the accelerations, but if the acceleration measurements are performed on the moving craft then no it's not possible alone by accelerometer readings.

Eventhough it's sufficient to use those accelerations $a_x$, $a_y$, and $a_z$ on a cartesian coordinate system to compute the displacement associated with the center of mass of the craft that's moving in 3D, nevertheless you will also need to know the orientation of the craft (hence the accelerometer) so that you can transform its acc measurement vectors into the navigation frame correctly.

Therefore, you also need to consider (at least) the gyroscopic readings to determine the craft orientation wrt your navigation frame. Once you determine the orienation of the craft, then you convert its current acceleration measurements into the non-rotating intertial frame in which you compute the track of the craft by integrating the accelerometer readings.


Understanding values saved in WAVE file The values that one can find in a wave file, e.g.

    0.0036
    0.0026
    0.0174
    0.0050
    0.0026
    0.0108
    0.0154
   -0.0114
   -0.0006
    0.0115

are - as I assume - relative values between -1 and 1, the range you normally see when watching a waveform in some audio editing program. Now, how are these relative values mapped when you playback a file? E.g. if there is a value of 0.5, does it mean playback at half of the maximal possible loudness? Is there some kind of reference value stored in the WAVE file to which the relative values refer to?

Another thing: In the example above there are 7 positive values follow by 2 negative values. So the membran vibrating in a speaker to create the sound is pressed longer in one direction than in the other direction, or how does one has to interpret this? Positive values are one direction of the vibration and negative ones are vibrations into another direction or am I wrong about this?

Thank you very much :-)

What you see are time-domain samples that correspond to the displacement of a speaker from some rest position as a function of time. For example, the first sample corresponds to the membrane's position at $t=0$, the next represents its position at $t=T$, and so on, assuming that there are $T$ seconds between samples. Common sampling intervals for audio include $T = \frac{1}{8000}$ sec and $T = \frac{1}{44100}$ sec.

In order to turn a stream of samples into real audio, a digital-to-analog converter (DAC) is used first to generate a time-varying analog voltage that corresponds to the samples. This analog voltage would be applied to the input of a power amplifier that drives the speaker. The speaker is designed in such a way that the position of its membrane varies proportionally to the speaker's input voltage. At the end of this chain, then, you have a speaker membrane that moves in accordance with the original audio samples, producing the desired audio via its vibrations.

As you noted, the sample values are normalized to a range of $\pm1$. The resulting loudness of the audio can be affected by many non-digital factors, such as the gain of the amplifier driving the speaker or the DAC's output voltage range. The peak sample values $\pm1$ just correspond to "maximum membrane displacement that the hardware can generate." Therefore, there is no way to encode absolute loudness of the audio in such a stream of samples, hence the use of such a normalization scheme.


Is there an optimal way to quantize log-likelihood ratios of an AWGN channel? My understanding is that given an AWGN channel and BPSK modulation, an LDPC decoder that uses message passing takes as input log-likelihood ratios $L$ of the following form by Bayes' rule:
$$
L=\frac{P(a=1|r)}{P(a=-1|r)}=\frac{2}{\sigma^2}r
$$
where $a\in\{\pm1\}$, $n\sim \mathcal{N}(0,\sigma^2)$, and $r=a+n$.
I have been referencing page 22 of this report, which uses a uniform quantizer and mentions that a larger dynamic range results in larger step sizes assuming a fixed number of bits.
Are there other quantization methods that might be better suited to the Gaussian nature of this channel? Said differently, are there quantization methods that better account for $\sigma^2$?

It depends on your actual decoder architecture, but if memory serves me right, uniform quantization is optimal for Min-Sum or Sum-Product decoders. (Intuitively, kind of makes sense – any other quantization could be modelled as "let's first apply a function to the LLRs and then quantize the result of that uniformly", and LLRs are what asymptotically achieves maximum likelihood performance.)
The question is where you put your quantization boundaries: If you have a fixed number of bits for your LLRs, do you quantize the region from, say, $-2\sigma$ to $+2\sigma$, because everything larger is sufficiently "safely" the sign-indicated value, anyway? Or do you have a large code word length, and you want to make sure that in sea of large LLRs, where some represent bit errors, the "righter" ones are
still understood as such?
That's a question of looking your design SNR and your target FER; no general answer or formula is known, far as I'm aware.


$\int_{-\infty}^{+\infty} |G(f)| \,e^{j2\pi ft}df=|g(t)|$? Given the absolute value of the Fourier transform of a signal $g(t)$:
$|G(f)|$
If I compute the inverse Fourier transform of $|G(f)|$,
$$\int_{-\infty}^{+\infty} |G(f)|\,  e^{j2\pi ft}df$$
do I obtain $|g(t)|$?

The answer is clearly 'no', except for trivial cases. Note that
$$G(f)=|G(f)|e^{j\phi(f)}\tag{1}$$
where $\phi(f)$ is the phase of $G(f)$. Rewriting $(1)$ as
$$|G(f)|=G(f)e^{-j\phi(f)}\tag{2}$$
and taking the inverse Fourier transform gives
$$\tilde{g}(t)=(g\star a)(t)\tag{3}$$
where $\tilde{g}(t)$ is the inverse Fourier transform of $|G(f)|$, and $a(t)$ is the impulse response of a (generally non-causal) all-pass with frequency response $A(f)=e^{-j\phi(f)}$.
If we are given $|G(f)|$, we do not know $\phi(f)$ and, consequently, we don't know $a(t)$, so we can't compute $\tilde{g}(t)$ from $(3)$. Anyway, $(3)$ shows the relation between the inverse Fourier transform of $|G(f)|$ and $g(t)$.
Note that since $|G(f)|$ is real-valued, $\tilde{g}(t)$ must be even
$$\tilde{g}(t)=\tilde{g}(-t)\tag{4}$$