Energy normalization across wavelet subbands According to the following reference: A Really Friendly Guide to Wavelets, © C. Valens, 1999
Equation 3, the wavelets are generated from the mother wavelet by scaling and translation. S is the scaling factor. The equation is divided by the square root of S for energy normalization across the different bands. I can not understand why this division achieves energy normalization.

According to the same reference, equation 14,  From Fourier theory we know that compression in time is equivalent to stretching the spectrum and shifting it upwards. I can understand spectrum stretching in this equation, but where is the spectrum shifting in the equation?1

The report under reference is A Really Friendly Guide to Wavelets by Clemens Valens, 1999.
The title of the question might be misleading. The energy normalization concerns the wavelets themselves, required to have the same energy, whatever the scaling and the shift. This is useful for theoretical and practical reasons (algorithmic/computational stability,  and several statistical properties).
If function $\psi(t)$ is of finite energy $\|\psi\|_2^2$ (or square integrable), then by a change of variable $t\mapsto st$, you will find that  $\psi\left(\frac{t-\tau}{s}\right)$ has finite energy as well, and that $\|\psi\left(\frac{t-\tau}{s}\right)\|_2^2 = |s|\|\psi\|_2^2$. For the wavelets to have the same energy, you should divide by $1/\sqrt{|s|}$.
An  interpretation is the following. The energy of function is homogeneous to a rectangle with a time span horizontally, and a square in amplitude vertically.
You want to keep the area of the rectangles constant for each wavelet. So if you dilate the time by $s>0$, you should contract the square in amplitude by $1/s$, then the amplitude by $1/\sqrt{s}$.
For question 2, the spectrum will be stretched/shifted accordingly, if you compute the Fourier transform of $\psi\left(\frac{t-\tau}{s}\right)$. The scaling property of Fourier is mentioned in Clemens Valens tutorial at Equation 14:
$$\mathcal{F}\left(f(at)\right) = \frac{1}{|a|}F\left(\frac{\omega}{a}\right)\,.$$


ICA - Statistical Independence & Eigenvalues of Covariance Matrix I am currently creating different signals using Matlab, mixing them by multiplying them by a mixing matrix A, and then trying to get back the original signals using FastICA.

So far, the recovered signals are really bad when compared to the original ones, which was not what I expected.

I'm trying to see whether I'm doing anything wrong. The signals I'm generating are the following:

s1 = (-x.^2 + 100*x + 500) / 3000; % quadratic
s2 = exp(-x / 10); % -ve exponential
s3 = (sin(x)+ 1) * 0.5; % sine
s4 = 0.5 + 0.1 * randn(size(x, 2), 1); % gaussian
s5 = (sawtooth(x, 0.75)+ 1) * 0.5; % sawtooth

One condition for ICA to be successful is that at most one signal is Gaussian, and I've observed this in my signal generation.

However, another condition is that all signals are statistically independent.

All I know is that this means that, given two signals A & B, knowing one signal does not give any information with regards to the other, i.e.: P(A|B) = P(A) where P is the probability.

Now my question is this: Are my signals statistically independent? Is there any way I can determine this? Perhaps some property that must be observed?

Another thing I've noticed is that when I calculate the eigenvalues of the covariance matrix (calculated for the matrix containing the mixed signals), the eigenspectrum seems to show that there is only one (main) principal component. What does this really mean? Shouldn't there be 5, since I have 5 (supposedly) independent signals?

For example, when using the following mixing matrix:

A =

0.2000    0.4267    0.2133    0.1067    0.0533
0.2909    0.2000    0.2909    0.1455    0.0727
0.1333    0.2667    0.2000    0.2667    0.1333
0.0727    0.1455    0.2909    0.2000    0.2909
0.0533    0.1067    0.2133    0.4267    0.2000

The eigenvalues are: 0.0000 0.0005 0.0022 0.0042 0.0345 (only 4!)

When using the identity matrix as the mixing matrix (i.e. the mixed signals are the same as the original ones), the eigenspectrum is: 0.0103 0.0199 0.0330 0.0811 0.1762. There still is one value much larger than the rest..

Thank you for your help.

I apologise if the answers to my questions are painfully obvious, but I'm really new to statistics, ICA and Matlab. Thanks again.

EDIT

I have 500 samples of each signal, in the range [0.2, 100], in steps of 0.2, i.e. x = 0:0.1:100.

Also, given the ICA Model: X = As + n (I'm not adding any noise at the moment), I am referring to the eigenspectrum of the transpose of X, i.e. eig(cov(X')).

UPDATE

As suggested (refer to comments), I tried FastICA on only 2 signals. The results were quite good (see pic below). The mixing matrix used was A = [0.75 0.25; 0.25 0.75]. However, the eigenspectrum  0.1657 0.7732 still showed only one main principal component.

My question therefore boils down to the following: What function/equation/property can I use to check whether a number of signal vectors are statistically independent?

Signals 3 and 5 appear to be quite correlated -- they share their first harmonic. If I were given two mixtures of those, I wouldn't be able to separate those, I'd be tempted to put the common harmonic as one signal and the higher harmonics as a second signal. And I would be wrong! This might explain the missing eigenvalue.

Signals 1 and 2 do not look independent too.

A quick and dirty "sanity-check" for the independence of two series is to do a (x, y) plot of one signal against the other:

plot (sig3, sig5)

and then to do the same (x, y) plot with one signal shuffled:

indices = randperm(length(sig3))
plot(sig3(indices), sig5)

If the two plots have a dissimilar look, your signals are not independent. More generally, if the (x, y) plot of the data shows "features", disymmetries etc, it's a bad omen.

Proper tests for independence (and those are the objective functions used in the ICA optimization loop) include, for example, mutual information.

ICA is recovering the most independent signals a linear mixing of which yields your input data. It will work as a signal separation methods and recover the original signals only if those were maximally independent according to the optimization criterion used in your ICA implementation.


how to convert SPL in dBA to dB I wonder how I can convert dBA (A weighted) to dB and vice versa.
can anybody provide any kind of table or online conversion app?
thanks
following formula comes from https://en.wikipedia.org/wiki/A-weighting
:\begin{align}
  R_A(f) &= {12194^2 f^4 \over \left(f^2 + 20.6^2\right)\ \sqrt{\left(f^2 + 107.7^2\right)\left(f^2 + 737.9^2\right)}\ \left(f^2 + 12194^2\right)}\ ,\\[3pt]
    A(f) &= 20\log_{10}\left(R_A(f)\right) - 20\log_{10}\left(R_A(1000)\right) \\
         &\approx 20\log_{10}\left(R_A(f)\right) + 2.00
\end{align}
how can I inverse this formula?
I also must mention that I have a specific frequency and I just want to convert this specific frequency SPL

You can't convert these directly. You need to start with either a calibrated pressure time domain pressure wave form or a calibrated power spectrum. To create calculate unweighted dBSPL yo simply sum the total energy. To calculate dBA, you need to apply an A-weighting filter to either one of those first.


How do I test stability of a MIMO system? Let's say I have a system similar to two interconnected IIR filters described like this:

\begin{align}
x_1(t)&=a_{11} x_1(t-1)+a_{12} x_1(t-2) +a_{13} x_2(t-1) + a_{14} x_2(t-2)+y_1(t)\\
x_2(t)&=a_{21} x_2(t-1)+a_{22} x_2(t-2) +a_{23} x_1(t-1) + a_{24} x_1(t-2)+y_2(t)
\end{align}

where $x$ are the outputs and $y$ are the inputs.

How would I go about testing the stability of such a system?

I myself have pretty much no background in MIMO systems, so I'm not really sure where even to start.

For a SISO system, I'd start with poles and zeros and go on from there, but I'm not sure how such a graph might looks like in the MIMO case.

Let's try to make your system into a state space representation with state:
$$
\mathbf{x}(t) = \left [
\begin{array}{c}
x_1(t)\\
x_2(t)\\
x_1(t-1)\\
x_2(t-1)\\
\end{array}
\right ]
$$
so that
$$
\mathbf{x}(t) = \mathbf{A} \mathbf{x}(t-1) + \mathbf{y}(t)
$$
where
$$
\mathbf{A} =
\left [ \begin{array}{cccc}
a_{11} & a_{12} & a_{13} & a_{14}\\
a_{21} & a_{22} & a_{23} & a_{24}\\
1   & 0 & 0 & 0\\
0   & 1 & 0 & 0\\
\end{array}
\right ]
$$
and
$$
\mathbf{y}(t) = \left [
\begin{array}{c}
y_1(t)\\
y_2(t)\\
0\\
0
\end{array}
\right ]
$$
or you could make $\mathbf{y}$ $2 \times 1$ and introduce a $\mathbf{B}$ matrix.

The system may be unstable if the eigenvalues of $\mathbf{A}$ are outside the unit circle.


How do I implement an adaptive thresholding algorithm for underwater sonar I want to implement an adaptive thresholding algorithm in MATLAB for filtering data received by an underwater sonar receiver. The data received has an interactive noise component resulting from underwater noise and specular reflection. CFARD method is close, but it does not serve my purpose. I have to image the data so that i am able to see the object, on a screen, that is placed underwater within the scandepth of the sonar. Any help will be greatly appreciated.

EDIT:

It is an underwater environment. I am trying to threshold a signal that has been received from a sonar transducer after it has been reflected by a solid target, that is situated in the same environment as the transducer. The problems belongs to Underwater Acoustic Imaging sonar domain. Problem is, I haven't been able to model the underwater environmental noise. From what I have read till now about this topic, the noise model follows a $K$-distribution. Also the environmental noise is not additive in nature, rather it's interactive. Hence the threshold has to be adaptive. I also mentioned CFARD method in my question. That is useful for signal processing in radar applications as we're just interested in finding a single point in a large area which has high energy. Same cannot be said about underwater acoustic imaging sonar, where we try to display the target on the screen as a video. I hope I have made it more clear now.

Your question has received quite few contributions, probably because of a lacking content. During a recent conference , I came across the PhD thesis: Détection en Environnement non Gaussien (Detection in a non-Gaussian environment). Since it is in French, I reproduce the abstract here:

  For a long time, radar echoes coming from the various returns of the
  transmitted signal on many objects of the environment (clutter) have
  been exclusively modelled by Gaussian vectors. The related optimal
  detection procedure was then performed by the classical matched
  filter. Then, the technological improvement of radar systems showed
  that the true nature of the clutter could not be considered as
  Gaussian any more. Though the optimality of the matched filter is no
  more valid in such cases, CFAR techniques (Constant False Alarm Rate)
  were proposed for this detector in order to adapt the value of the
  detection threshold to the multiple local variations of the clutter.
  In spite of their diversity, none of these techniques turned to be
  either robust or optimal in these situations. With the modelling of
  the clutter by non-Gaussian complex processes, such as SIRP
  (Spherically Invariant Random Process), optimal structures of coherent
  detection have been found. These models describe many non-Gaussian
  laws, like K-distribution or Weibull laws, and are acknowledged in the
  literature to model many experimental situations in a relevant way. To
  identify the law of their characteristic component (namely the
  texture) without statistical a priori on the model, we propose, in
  this thesis, to tackle the problem by a Bayesian approach. Two new
  estimation methods of the texture law emerge from this proposition:
  the first one is a parametric method, based on a Padé approximation of
  the moment generating function, and the second one results from a
  Monte Carlo estimation. These estimations are carried out on reference
  clutter data and lead to two new optimal detection strategies,
  respectively named PEOD (Padé Estimated Optimum Detector) and BORD
  (Bayesian Optimum Detector Radar). The asymptotic expression of the
  BORD (convergence in law), called the "Asymptotic BORD", is
  established together with its law. This last result gives access to
  the optimal theoretical performances of the Asymptotic BORD, and may
  also applied to the BORD if the data correlation matrix is
  non-singular. The detection performances of BORD and those of
  Asymptotic BORD are evaluated on experimental ground clutter data. We
  obtained results that validate both the relevance of SIRP model for
  the clutter, the optimality of the BORD and its adaptability to any
  type of environment.

The math should be readable. If it is of any help, you can track English reference by the author or PhD thesis committee.


How to design bandpass filter in python when centre frequency is greater than sampling frequency? I implemented basic bandpass filter using Scipy resource. But how to implement the same when signal is centred at non-zero centre frequency with centre frequency > sampling frequency.

I trying to implement bandpass filter as current code throws error

  ValueError: Digital filter critical frequencies must be $0 \le \omega_n \le 1$

which clearly makes sense.

Example:

Centre frequency: $100\textrm{ MHz}$
Sampling frequency: $5\textrm{ MHz}$
Highcut: $102\textrm{ MHz}$

Clearly (highcut/nyq) > 1

As far as the digital signal is concerned, the signal only exists uniquely in one Nyquist interval (which is $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate. Every other interval, such as in your range of interest that is well above the sampling rate, is a replica of the primary interval.  Given that, note that the range from $f_s/2$ to $f_s$ is the same as $-f_s/2$ to $0$, so we can also describe the primary Nyquist interval as being from $0$ to $f_s$ (and Matlab's fftshift command converts a frequency axis from $0$ to $f_s$ to $-f_s/2$ to $+f_s/2$ according to this relationship). The following picture helps illustrate this:

This is formed because sampling is multiplying the time domain signal by an impulse train. An impulse train in time is an impulse train in frequency (as represented by the second row in the figure above). Multiplying in time is convolving in frequency, so the result is to convolve our analog spectrum with the spectrum of the sampling process (an impulse train) resulting in the replicas of the spectrum in the Digital Output Spectrum as shown. Since digitally every higher frequency is an exact replica of what is given from $-f_s/2$ to $+f_s/2$, we need only concern ourselves with the frequency over that range, which is the unique freqency span of a digital signal.  However what I refer to as the "unwrapped view" of the digital signal has helped give me a much more intuitive insight into many multirate signal processing challenges, as well as undersampling (or bandpass sampling) as discussed here. The "unwrapped view" is done by extending the digital frequency axis from $-\infty$ to $+\infty$ while maintaining mathematical equivalence, which means replicating the same signal in each unique span of $n[0:f_s]$ where n is an integer extending to $-\infty$ to $+\infty$ as I have done in the figures here for the "Digital Output Spectrum".

That said, you would need to map the signal of interest to the primary interval and then design the filter in that space.

Further if your signal is real, the positive and negative frequencies would be complex symmetric (as diagrammed also in my figure), in which case we can just concern ourselves with the signal from $0$ to $f_s/2$. I suspect but have not confirmed the range in the python tool you referenced, given as 0 to 1 refers to DC to $f_s/2$ specifically.

As you suggest, we can do bandpass sampling as long as the A/D Converter has an analog input bandwidth that supports this, and I won't get into the detail of it here except to say that the jitter requirements on your sampling clock are more stringent with bandpass sampling. This implies you have an appropriate bandpass filter (refered to as the anti-alias filter) in your analog system prior to your A/D Converter, as the A/D Converter will downconvert every signal in every Nyquist zone if unfiltered energy exists (and within the analog bandwidth of the A/D Converter, which effectively becomes part of the anti-alias filter). This bandpass filtering operation is depicted in the figures below, showing how the energy in the other Nyquist zones gets translated to the same digital spectrum. Without analog bandpass filtering, all these bands will downconvert through the sampling process to the same frequency range, leading to a significant increase (increase in the noise floor even if there are no other signals present if the amplified noise floor in front of the A/D Converter is not bandpass or lowpass filtered in any way, as said earlier, out to the input analog bandwidth of the A/D Converter.

That said, in your case of a 100 MHz center frequency, the 100 MHz lands directly on a harmonic of the sampling clock and therefore maps to 0. Unless you are dealing with a complex signal or a real signal in a synchronous system, this would be a poor choice for a center frequency or sampling rate (by synchronous, I mean the signal of interest is coherent to the sampling clock such that you are assured it is exactly centered on 100 MHz and complex symmetric on each side of 100 MHz). The reason is, any offset in center frequency will cause aliasing in the image that gets mapped to 0 resulting in severe distortion. (Exceptions to this would be in closed-loop systems that will drive any offset to exactly 0). If you use complex sampling with resulting quadrature I and Q signal paths then the spectrum is unique from $-f_s/2$ to $+f_s/2$ and this will be a non-issue. For real signals it is common practice to center the signal near $f_s/4$ and multiples there-of just as I have done in these figures as it optimally relaxes the requirement on your analog anti-alias filter (the images will be equidistant as in my figures). It is common to then use a complex down-converter to then move the signal to baseband appropriately without aliasing issues. So to do this you would either shift your center frequency to 102.5 MHz, or shift your sampling clock (using some math) such that 100 MHz is mid-way in a Nyquist zone. You would then design a bandpass digital fitler that is in your primary Nyquist zone of $0$ to $f_s/2$ or as your tool is givng you in normalized terms $0$ to $1$.


Numpy's real fft (rfft) - losing power Related to another problem I'm having, I was looking into the workings of numpy's rfft2 and irfft2. These are special versions of the FFT routine, in so far that it needs less input; because you require the real-space image to be real you only need to 'fill' half of Fourier space - due to symmetry, that's all the information you need.

Now, in the docs it is mentioned that performing the inverse transform right after the forward transform should return the original array - to within numerical precision. Time to put this to the test.

If I do the following:

import numpy as np

fourier_image = np.random.uniform(-1,1,(100,51))
fourier_recovered = np.fft.rfft2(np.fft.irfft2(fourier_image))

print 'Original power: ', sum(sum(abs(fourier_image)**2.))
print 'Recovered power: ', sum(sum(abs(fourier_recovered)**2.))

... I find the original power and the recovered power differ at the percent level. That seems a rather big difference to me. Am I doing something wrong, or is this just the best we can do?

I naively assumed any (N,N/2+1) Fourier image had the right symmetry, but this goes wrong at the 0- and Nyquist frequency because they don't have a symmetric counterpart. Here, you have to impose the symmetry explicitly.

So, if we generate an image:

N = 100
fourier_image = np.random.uniform(-1,1,(N,1+N/2))

we impose the required symmetry as follows:

fourier_image[1:N/2,0] = np.conj(fourier_image[N:N/2:-1,0])
fourier_image[1:N/2,-1] = np.conj(fourier_image[N:N/2:-1,-1])

And doing the same as in the question above, we find the input - and recovered power:

print 'Original power: ', sum(sum(abs(fourier_image)**2.))
print 'Recovered power: ', sum(sum(abs(fourier_recovered)**2.))

which are indeed the same:

Original power: 1664.3177945
Recovered power: 1664.3177945

For an odd-shaped image this will be slightly different but the idea is the same.


What do colored noises look like in the time domain? I understand that the time domain representation of white noise looks like impulses. How do colored noises like brown, pink etc. look like when we perform an inverse Fourier Transform on them ? What could be some sources of colored noise that may affect speech signals?

White noise implies no correlation between samples of the noise, even consecutive samples.  Colored noise, therefore, implies that there is correlation of some sort between the noise samples, which in turn implies that we can take advantage of that correlation to get rid of some of the noise.

Beyond that, there is not a lot that we can say about what it looks like in the time domain other than to say that colors that favor low frequencies, like pink noise, will tend to change more slowly, while high frequency colors, like "blue" noise, will tend to change more quickly.


Inverse Problem / Deconvolution with Pink Noise Hi I dived somewhat into deconvolution of systems which can be described as:

$s(t) = o(t) * h(t) + n(t)$

where $s$ is my measured 1D time resolved signal, $o$ is the original signal $h$ is the kernel or PSF pf my system, which is convolved with the original signal (with $*$ beeing the convolution operator), $n$ as sample-function for $1/f^\alpha$ noise ($\alpha > 0$).

For all I know there are several methods to deconvolve for an estimate of o. When choosing the correct method, the consideration of the noise is important. E.g. I should use the wiener deconvolution for systems with white noise, lucy-richardson for possion noise.

Which method should I use for pink noise (given the fact that there are several other (read: deeper mathematical) aspect which should be taken into account)?

On a side note: any recomendations for Papers/Books on this topic (Deconvolution of systems as described)? Thanks in advance!

Edit:
For the system I know an estimate of the kernel function $h$ and I can make blank meassurements to estimate the nature of the nois sample function (it ssemst to be white noise for small freqencies and pink noise for higher freqencies).

Some Remarks:

2nd Order Knowledge
It seems you only have a knowledge about the Auto Correlation of your data (2nd Order) and not on its distribution. Hence methods you can apply are ones which minimizes only 2nd order functions of the noise (Such as MMSE).
Wiener Filter Minimizes the MMSE
As can be seen in the derivation of the Wiener Deconvolution Filter the only assumption made is independence from the signal. If that holds in your data this is a good solution. Pay attention that Wiener Filter model is that the input data is also stochastic and all data is wide sense stationary.
Richardson Lucy Deconvolution
Pay attention this is a Maximum Likelihood method. Namely it assumes the data is a parameter, not a random data. It also assumes the distribution of the noise and not only its second moment,
Decide the Cost Function
First thing you should decide what would you like to minimize according to your model. Does the MMSE makes sense? Maybe a Least Squares is good enough? May you come with other objective. This is the first step.


Sum of Sine and Cosine with Random Phase as LTI System I have the following system:

Where $ {H}_{1} \left( f \right) = {H}_{2} \left( f \right) $ and $ \theta \sim U[0, 2\pi]$ independent of any other factor in the system.

Given the input is identical, Is an LTI system?
Could you prove it?

Otherwise I don't have explanation why if the input is White Gaussian Noise with the same Variance to each of the inputs (Same statistics, not the same signal) $ {Y}_{1} $ and $ {Y}_{2} $ aren't Gaussian Process yet $ r(t) $ is a Gaussian Process.

I have the feeling the random phase is offset by the sum of the Sine and Cosine.

Thank You.

So let's just look at it from a linear systems perspective:
$$
\begin{array}{ccl}
r(t) &=& Y_1(t) + Y_2(t)\\
&=& Z_1(t) \cos(2\pi f_0 t + \theta) + Z_2(t) \sin(2\pi f_0 t + \theta)\\
&=& [ X_1(t) * H_1(t) ] \cos(2\pi f_0 t + \theta) + [ X_2(t) * H_1(t)] \sin(2\pi f_0 t + \theta)
\end{array}
$$
where $X_1$ and $X_2$ are the respective (unlabeled) inputs, $H_1$ and $H_2$ are the impulse responses of the filters, and $*$ is convolution.

You appear to be asking whether the system is LTI when $H_1 = H_2 = H$ and $X_1 = X_2 = X$.  In that case, we get
$$
\begin{array}{ccl}
r(t) &=& [ X(t) * H(t) ] [\cos(2\pi f_0 t + \theta) + \sin(2\pi f_0 t + \theta)]\\
&=& \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)
\end{array}
$$

So it's linear:
$$
\begin{array}{ccl}
r_a(t) &=& \sqrt{2} [ a(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_b(t) &=& \sqrt{2} [ b(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_{a+b}(t) &=& \sqrt{2} [ (a(t) + b(t)) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4) = r_a(t) + r_b(t)\\
\end{array}
$$
But it's not time-invariant:
$$
\begin{array}{ccl}
r(t) &=& \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_\delta(t) &=& \sqrt{2} [ X(t-\delta) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r(t-\delta) &=& \sqrt{2} [ X(t-\delta) * H(t) ] \cos(2\pi f_0 (t-\delta) + \theta - \pi/4)\\
\end{array}
$$
as $r_\delta(t) \not = r(t-\delta)$.

EDIT

So, what happens if $X_1 = X_2 = X$ and $X\ \tilde{\ }\ N(0,\sigma^2)$ and $\theta\ \tilde{\ }\ U[0,2\pi)$?

I can see why $Y_1(t) = [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta)$ might be non-Gaussian.

However, I cannot see why $Y_1 + Y_2 = \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)$ will be Gaussian?

Can you please update your question to explain?


About how to increase the FIR filter sampling frequency in FPGA? And what is the tradeoff of increasing the sampling frequency? I am trying to implement a FIR filter on FPGA and trying to have a solid understanding of the FIR filter tap delay and sampling frequency.

Does the “one tap” delay equal to “1/Fs (sampling frequency)”? If I have N-tap, the total delay will be N/Fs? If the Fs sampling frequency is increased, the “one tap” delay is decreased?

Is there any trade off/drawback of increasing the sampling frequency? Maybe it will take more processing time for the output y(n) to come out?

Thanks

Does the “one tap” delay equal to “1/Fs (sampling frequency)”?

Yes

  If I have N-tap, the total delay will be N/Fs?

Depends on how you define "total" delay, but in general the answer is no. For a minimum phase filter the delay will be 1. For a linear phase filter it would roughly be N/2

  If the Fs sampling frequency is increased, the “one tap” delay is decreased?

Yes. But it also changes the frequency response of your filter.

  Is there any trade off/drawback of increasing the sampling frequency?

Many. It really depends on what the requirements of your application are.


How to obtain the Sound Pressure Level from .wav file based on these given parameters? I have a wav file recorded from a hydrophone. I am not quite familiar with signal processing so I really need some advice from you.

I have gathered the information which looks useful and related.

sample rate = 250 000 Hz

Receive sensitivity of hydrophone = -180 db re 1V/μ Pa

Peak-to-peak voltage range = 5.64 V

Pre-amplifier gain = 18 dB (Modified)

Let A be the amplitude saved in the wav file. What is the formula to transform A from amplitude to Sound Pressure Level?

If the given information is not sufficient, please tell me what more do I need?

Any helps are appreciated!

The hydrophone has -180 dB 1V/uPa or -60 dB V/Pa. 1 Pa produces hence 1mV or, other way around, the sensitivity is 1mV/Pa. Peak voltage is 5.64 peak to peak. That's 2.82V peak or 2V RMS.

A full scale sine wave has 2 V RMS which is 2000 Pa. dBSPL is referenced against 20uPa so 2000 Pa is equivalent to 160 dBSPL.

The clipping point of your wave file corresponds to 160 dB SPL. The exact relationship between number and dBSPL  depends  on the format of the wave file (fixed vs. floating, number of bits). For a 16 bit fixed point format, the clipping point corresponds to the number 32768.


Looking for a simple and effective notch filter I ran a finite-difference simulation and the behavior of an output signal, $s$, in time, $t$ (sampled with period $\Delta t$) behaves approximately as in the figure below. It is well-described by a slowly-varying pulse plus a rapidly oscillating signal of the same envelope with lower amplitude.

For concreteness, I computed this sample signal in MATLAB with the following code.

f_noise = 20;
t = linspace(-4,4,2e3);
signal = exp(-t.^2);
noise = 0.25*exp(-t.^2).*sin(2*pi*f_noise*t);
signal_w_noise = signal + noise;

I emphasize that this is just a sample signal and not the actual data I want to analyze.

I wish to compute the "time-average" of the signal. By this I mean the signal without the rapidly oscillating component (which we could call noise). To get a clear picture of the problem, I calculated the Fourier transform of the signal approximately using MATLAB's fft function and obtained the spectra below.

From this plot, I suspect a notch-filter could do the trick. However, I know very little of signal processing and I know how to use first-order filters from undergraduate electric-circuit courses. What is a simple to use but effective notch-filter I could use? I am particularly looking for the transfer function in the frequency domain of such notch-filter, because the software I use to generate the data is equipped with a FFT algorithm, so I'd like to do the processing there and avoid exporting the whole data to, say, MATLAB.

i might also recommend a low-pass filter over a notch filter unless the OP is absolutely certain, in advance, what the frequency of the "noise" is.  if the OP has less information about the high-frequency noise that the OP desires to remove (leaving the moving average), then a LPF with DC gain of 1 (or 0 dB) and a cutoff frequency lower than the lowest expected frequency component in the noise is indicated.  not a notch filter.

whether it's a Kaiser-windowed-sinc FIR or an IIR LPF (perhaps a simple 2nd-order IIR suffices) is a design and performance issue.

but i would not recommend a notch unless it is known apriori that the "noise" is sinusoidal and what the frequency of the sinusoidal "noise" is.


How to Match 2 Signals which Are a Shifted and Scaled Version of Each Other I have 2 signals S1 and S2 that contain the same information, but S2 is shifted and scaled compared to S1 by an unknown amount (but small; eg shift would be of the order of 1-10 samples).

What is the best (=most effective yet simple) way to calculate (or estimate) the shift and scaling factor to apply to S2 in order to 'match' with S1.

I have read a few papers referring to using the phase information in FFT, but those were working on signals where scaling was not really an issue (and i must admit not to have understood them completely anyway). I am a bit lost on how to put this in practice when developing that in Matlab (actually I use Octave).

I would love to have an explanation using Matlab/Octave code if possible.

Thanks for any help.

Steve

PS: I add the tag 'image processing' because I suspect that such techniques would be used in automatic stitching of images (in panoramas for example)

Matlab has a function called xcorr() that computes the cross correlation between two vectors.  The cross correlation will be highest when the two signals overlap.  The following code will compute the shift in units of samples:

x = xcorr(s1,s2)
[val,idx] = max(x);
shift = length(s1) - idx

The scaling factor should be given by:

scale = max(s1)/max(s2)


Apply band pass filter (BPF) before calculating the power of a signal (PSD)? From my understanding PSD tries to give a good estimation on the power each frequency attributes to the overall signal power.

If I am only interested in a frequency range $(\omega_0 \leq \omega \leq \omega_1)$, could I put the signal through a band pass filter for that same frequency range then calculate the signal power using this formula:

$$P_x = \frac{1}{N_1 - N_0 + 1}\sum_{n = N_0}^{n=N_1} \left|x(n)\right|^2  $$

Would this technique give me what I am looking for, (the power of the signal from $N_0$ to $N_1$ for the specific frequency range)?

How does this compare to PSD estimation using any of the popular techniques?

Pros, Cons?

Thank you!

You can switch to the frequency domain with regards to Parseval's theorem. I guess you can transform your problem via discrete-time fourier transformation.


Guard Interval in OFDM If I extend the length of the Cyclic Prefix (Guard Interval) in an OFDM signal, the loss in data rate is clear by means of the time domain, but I can't see it in the frequency domain as the BW is given by Nsc/T_useful.

What am I missing?

I think that maybe calculating BW out of Nsc/T_useful is only correct after removing the CP in the demodulation process while the transmitted signal won't look the same "on the air", is that true?

There isn't any obvious effect in the frequency domain. However, extending the duration of the cyclic prefix does make your signal more resistant to multipath channels with a large delay spread. As a rule of thumb, you need to size the cyclic prefix duration to cover the entire delay spread. If this condition is met, then the linear convolution that is imposed on your signal by the multipath channel will manifest itself as circular convolution on the cyclic-prefixed signal. This circular convolution can be equalized by simple pointwise multiplication in the frequency domain at the receiver.

This is why one of OFDM's advantages is often cited as the fact that it permits the use of a "one-tap equalizer." In single-carrier systems, equalization can be one of the harder things to get right, so anything that can simplify that portion of the receiver is a help.


Phase plot of $H(e^{j\omega})=(1-re^{-j\omega})\left(1-\dfrac{1}{r}e^{-j\omega}\right)$, $0<r<1$ I want to find the phase plot of $H(e^{j\omega})=(1-re^{-j\omega})\left(1-\dfrac{1}{r}e^{-j\omega}\right)$, $0<r<1$ for the interval $0\leq \omega \leq \pi$.

Method 1:

$H(e^{j\omega})=1-\left(r+\dfrac{1}{r}\right)e^{-j\omega}+e^{-2j\omega}=e^{-j\omega}\left[2\cos \omega-\left(r+\dfrac{1}{r}\right)\right]$

$\angle H(e^{j\omega})=-\omega+\angle \left[2\cos \omega-(r+\dfrac{1}{r})\right]=-\omega+\pi$

Method 2:

$\angle H(e^{j\omega})=\angle(1-re^{-j\omega})+\angle(1-\dfrac{1}{r}e^{-j\omega})$

$$\angle H(e^{j\omega})=\begin{cases}
\tan^{-1}\bigg(\dfrac{r\sin \omega}{1-r\cos\omega}\bigg)+\pi-\tan^{-1}\bigg(\dfrac{\sin \omega}{\cos \omega-r}\bigg),0\leq \omega\leq \cos^{-1}(r)\\
\tan^{-1}\bigg(\dfrac{r\sin\omega}{1-r\cos\omega}\bigg)+\tan^{-1}\bigg(\dfrac{\sin\omega}{r-\cos \omega}\bigg),\cos^{-1}(r)\leq \omega\leq \pi\\
\end{cases}$$

I tried applying the formula $\tan^{-1}x+\tan^{-1}y=\tan^{-1}\bigg(\dfrac{x+y}{1-xy}\bigg)$, but I am confused when to use this formula and unable to proceed forward. Please help in solving this.

In fact you've already proved that the two expressions must be identical. But if you want to show it the hard way, you can indeed use the formula you suggested. The argument of the combined arctangent function becomes

$$\begin{align}\frac{\frac{r\sin(\omega)}{1-r\cos(\omega)}+\frac{\sin(\omega)}{r-\cos(\omega)}}{1-\frac{rsin^2(\omega)}{(1-r\cos(\omega))(r-\cos(\omega))}}&=\frac{r\sin(\omega)(r-\cos(\omega))+\sin(\omega)(1-r\cos(\omega))}{(1-r\cos(\omega))(r-\cos(\omega))-r\sin^2(\omega)}\\&=\frac{\sin(\omega)(1-2r\cos(\omega)+r^2)}{-\cos(\omega)(1-r\cos(\omega)+r^2)+r(\underbrace{1-\sin^2(\omega)}_{\cos^2(\omega)})}\\&=\frac{\sin(\omega)(1-2r\cos(\omega)+r^2)}{-\cos(\omega)(1-2r\cos(\omega)+r^2)}=-\frac{\sin(\omega)}{\cos(\omega)}=-\tan(\omega)\end{align}$$

Now you just have to be careful with the principal value and you'll obtain the same result as with your first method (which is of course much smarter).


Hilbert transform from analytic signal Show that the Hilbert transform of $h(t) = m(t) \cos(2 \pi \nu_c t)$ is

$$\hat{h} (t) = m(t) \sin(2 \pi \nu_c t),$$

where $m(t)$ is a real valued, band-limited function (i.e. we have Fourier transform $M(\nu) = 0$ for $|\nu| > \nu_m$) and $\nu_c > \nu_m.$

Attempt

I did some research and I found out that this result follows directly from 'Bedrosian's theorem'. But I am required to compute this by first finding its analytic signal $h_a$ whose imaginary part would then be the Hilbert transform. Here is my expression for $h_a (t):$

$$h_a (t) = 2 \int^\infty_0 \Big[ H(\nu) \Big] e^{j 2 \pi \nu t} \ d\nu = 2 \int^\infty_0 \Big[ \frac{1}{2} (M(\nu + \nu_c) + M(\nu - \nu_c)) \Big] e^{j 2 \pi \nu t} \ d\nu.$$

I have used the 'modulation property' of the Fourier transform to get to the RHS. So how can I proceed with the integration when we do not have an explicit expression for $m(t)$?

Note that the condition on the relation between the bandwidth of $m(t)$ and the modulation frequency ensures that $h(t)$ is a band pass signal. It can be written as

$$h(t)=\frac12 m(t)e^{j2\pi \nu_c t}+\frac12 m(t)e^{-j2\pi \nu_c t}\tag{1}$$

From our assumption $\nu_m<\nu_c$ it follows that the two expressions on the right-hand side of $(1)$ do not overlap in the frequency domain. Consequently, the analytic signal is simply given by the first term on the right-hand side of $(1)$ (scaled by a factor of $2$), which is the part of the spectrum of $h(t)$ occurring at positive frequencies:

$$h_a(t)=m(t)e^{j2\pi\nu_c t}\tag{2}$$

The Hilbert transform of $h(t)$ is given by the imaginary part of $(2)$:

$$\hat{h}(t)=\text{Im}\{h_a(t)\}=m(t)\sin(2\pi\nu_c t)\tag{3}$$

EDIT: The equation for $h_a(t)$ given in your question is also correct. If you want to continue that way, it's important to see that due to $\nu_c>\nu_m$, $M(\nu+\nu_c)$ equals zero for $\nu>0$. So the integral can be written as

$$h_a(t)=\int_0^{\infty}M(\nu-\nu_c)e^{j2\pi\nu t}d\nu\tag{4}$$

Likewise, $M(\nu-\nu_c)$ is zero for $\nu<0$, so we can extend the lower integration limit in $(4)$ to $-\infty$:

$$h_a(t)=\int_{-\infty}^{\infty}M(\nu-\nu_c)e^{j2\pi\nu t}d\nu=\mathcal{F}^{-1}\{M(\nu-\nu_c)\}=m(t)e^{j2\pi\nu_c t}\tag{5}$$

which is of course identical to $(2)$.


Frequency Domain Correction of Frequency Offset or Doppler Frequency in LTE We are working on LTE Uplink CFO Correction for 3 MHz band. The time domain CFO correction is easy, which depends on the processing of 3840 samples step by step multiply by $\exp(-j 2\pi f_\text{offset}t)$. . What will be the equivalent frequency domain correction?

Generally, a modulation with a complex sinusoid in time domain is equivalent to a shift in frequency domain (I find that quite intuitive, since it's the very definition of frequency shift).

That means that for $f_\text{offset} = M\frac{f_\text{sample}}{N_\text{FFT}},\quad M\in\mathbb Z$ you can just correct things by shifting the frequency domain vector left or right.

However, you also have to deal with the fractional frequency offset. If you don't, you'll get Inter-Carrier-Interference (ICI).

The only way to solve that is by interpolating the correct values from the values you have with sinc interpolation.

Sinc interpolation is easiest done by taking the IDFT, then modulating with a $e^{j2\pi f_\text{offset fractional} t}$, and then DFT'ing back.

In other words: that frequency correction is done in time domain.

So, do your fine frequency correction in time domain, not in frequency domain. Doing it in frequency domain requires at least the effort of going back to time domain.


Confused about shifting convolutions We know that for two signals $x[n]$ and $h[n]$ such that :
$$
y_{1}[n]=(x*h)[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]=\sum_{k=-\infty}^{\infty}h[k]x[n-k]
$$
We can deduce that :$$y_{2}[n]=x[n+2]*h[n]=y_{1}[n+2]$$ and $$y_{3}[n]=x[n]*h[n+2]=y_{1}[n+2]$$ but can we know what :
$$
y_{4}[n]=x[n+2]*h[n+2]
$$
gives in terms of $y_{1}[n]$? I am unaware if I might be asking a really stupid question or am I not seeing it? Is it $y_{1}[n+4]$?

Whenever you have doubts regarding the properties of the convolution operator, you should resort to its definition.
Let $x_s[n] = x[n+2]$ and $h_s[n] = h[n+2]$. Then:
$$y_4[n] = x_s[n] * h_s[n] = \sum_{k=-\infty}^\infty x_s[k]h_s[n-k]$$
If we replace with the definitions of our discrete functions:
$$y_4[n] = \sum_{k=-\infty}^\infty x[k+2]h[n-k+2]$$
See what happens if we make the change of variables $m = k+2$:
$$y_4[n] = \sum_{m=-\infty}^\infty x[m]h[n-(m-2)+2]=\sum_{m=-\infty}^\infty x[m]h[n+4-m]$$
Does that final expression sound familiar? You were right indeed in your OP, as that is exactly what you thought:
$$y_4[n] = \sum_{m=-\infty}^\infty x[m]h[n+4-m] = y_1[n+4]$$


Applications or physical interpretation of auto-convolution? I wonder if anyone has any experience with auto-convolution. In particular i'm interested in understanding the physical interpretation of it.
I understand what convolution, correlation and auto-correlation are, also i'm aware that the definition of auto-convolution will be something like
$$f\ast f = \int_{-\infty}^{\infty} f(\tau)f(t-\tau)d\tau $$
but i still don't get what are the implications o the meaning of it.  I've been looking for a while and so far i haven't found any good or detailed explanation (on constrast with auto-correlation).
So, if anyone has any experience dealing with this topic or has any intuitive interpretation that could share, i'll appreciate it. Thanks.

Autoconvolution is used in signal detection, but the way you've written it is not correct. Suppose you're trying to detect a signal $f(t)$ by filtering with h(t).

$y(t) = f(t) \ast h(t)$

You want to maximize your response to the signal $f(t)$. We can do this by maximizing the correlation coefficient between $f$ and $h$. Here the correlation is time-varying , so we'll maximize the average autocorrelation coefficient. We'll assume the $f$ and $h$ signals have DC values of zero for simplicity.I'll use $\mu_y$ to denote average value of a the auto-correlation of response, $y$.

$ h = argmax_{h} \ \   \mu_y\big(E_f E_h \big)^{-\frac{1}{2}}$

Let's take a differential of our performance $ J \propto \mu_y\big(E_f E_h \big)^{-\frac{1}{2}}$ with respect to h

$\partial_{h(\tau)} J \propto \partial_{h(\tau)} \bigg((\mu_y\big(E_f E_h \big)^{-\frac{1}{2}}\bigg)$

$\ \ \ \ \ \ \ \ \  =  \frac{1}{|dom(f)|}\bigg(\partial_{h(\tau)}\mu_y\bigg)\big(E_fE_h)^{-\frac{1}{2}} +\frac{1}{2}\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\bigg(\partial_{h(\tau)}E_h\bigg)$

$\ \ \ \ \ \ \ \ \  \propto  \bigg(\partial_{h(\tau)}\frac{1}{|dom(f)|}\int_{\infty}^{\infty}f(t-\tau)h(\tau)d\tau \bigg)\big(E_fE_h\big)^{-\frac{1}{2}} +\frac{1}{2}\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\bigg(\partial_{h(\tau)}E_h\bigg)$

$\ \ \ \ \ \ \ \ \  =  \bigg(\frac{1}{|dom(f)|}\int_{\infty}^{\infty} f(t-\tau)d\tau \bigg)\big(E_fE_h\big)^{-\frac{1}{2}} -\frac{1}{2}\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\bigg(\partial_{h(\tau)}E_h\bigg)$

$\ \ \ \ \ \ \ \ \  =  \big(E_fE_h\big)^{-\frac{1}{2}} \frac{1}{|dom(f)|}\int_{\infty}^{\infty} f(t-\tau)d\tau-\frac{1}{2}\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\bigg(\partial_{h(\tau)}\int_{\infty}^{\infty} h^2(\tau)d\tau \bigg)$

$\ \ \ \ \ \ \ \ \  =  \frac{1}{|dom(f)|}\big(E_fE_h\big)^{-\frac{1}{2}} \int_{\infty}^{\infty} f(t-\tau)d\tau-\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\int_{\infty}^{\infty} h(\tau)d\tau$

$\ \ \ \ \ \ \ \ \  \propto  \frac{1}{|dom(f)|}\int_{\infty}^{\infty} f(t-\tau)d\tau-\mu_yE_h^{-1}\int_{\infty}^{\infty} h(\tau)d\tau$

$\ \ \ \ \ \ \ \ \  = \int_{\infty}^{\infty} \bigg(\frac{1}{|dom(f)|}f(t-\tau)-h(\tau)\mu_yE_h^{-1}\bigg)d\tau$

At minimum $J$ we shouldn't assume $h(\tau) = 0$, we'll enforce $$\frac{1}{|dom(f)|}f(t-\tau)-h(\tau)\mu_yE_h^{-1} = 0$$

and pretty readily you get

$$h(\tau) = \frac{E_h}{\mu_y|dom(f)|}f(t-\tau)$$

Or most importantly

$$h(\tau) \propto f(t-\tau) \text{where} \ \ \tau \ \ \text{is time and} \ \ t \ \ \text{is the delay}$$

or, using perhaps better "variable names"

$$h(t) \propto f(t_d-t) \text{where} \ \ t \ \ \text{is time and} \ \ t_d \ \ \text{is the delay}$$

That is, if we want to maximize the correlation between our signal detector with impulse h(t), we better pick h(t) to be a time-reversed and time-shifted version of our signal of interest. In practice $t_d$ would probably be set to zero, as it just represents whenever the $f$ part that you're looking for finally arrives.

Under this chosen $h(t) \propto f(t_d - t)$, your original question makes more sense. The autocorrelation signal becomes proportional to the accumulated energy of $f(t)$ that is seen by your filter.

$y(t) = f(t) \ast h(t)$

$y(t) \propto \int f(\tau) h(t-\tau)d\tau$

$y(t) = \int f(\tau) f(t_d-(t+\tau))d\tau$

$y(t) = \int f(\tau) f(\tau + t_d - t)d\tau$

$y(t) = \int P_f(\tau + t_d - t)d\tau$


Edge detection in a crack pattern image Here's an image of a cracked pattern. I wish to convert it to binary image where the continuous crack lines form the only discernible pixels. The main problems I am having in the lack of uniform lightning, which confuses threshold tools in image processing software, and the "black dots" which should be filtered out.

I am trying to use either Gimp, imagemagick or ImageJ to accomplish this, but have so far not been able to get to decent result.

I have a partial solution, without filtering the black spots.

These are my suggestions. They are pretty low tech, so feel free to play around with them.
First get rid of the uneven illumination using a Bandpass Filter: These are the specific operations that I applied to your image recorded as a FIJI / ImageJ macro:

//convert to 32 bit image
//manually by: Image->Type->32 bit
run("32-bit");
//Apply bandpass Filter:
//Manually by: Process->FFT->Bandpass Filter...
run("Bandpass Filter...", "filter_large=75 filter_small=0 suppress=None tolerance=0");
//Enhance Contrast
run("Enhance Contrast", "saturated=0.35");

What that does is get rid of the slowly varying background in your image. Now you can manually adjust the threshold using Image->Adjust->Threshold. Press Apply after adjusting the threshold to your needs. This is what it looks like:

Again, this is all very low tech, but try to play around with this. Finally you might want to run Process->Binary->Skeletonize on the image.


How to derive Q from r in this filter definition? I am trying to build a resonance filter with real-time control of both center frequency and Q. I've come across what looks to me like a suitable "recipe" for a second-order IIR filter that might do the job, at http://www.music.mcgill.ca/~gary/307/week2/filters.html (toward the bottom of the page, quoted here with some loss of formatting):

$$y[n] = x[n] - a_1 y[n-1] - a_2 y[n-2]$$

where $a_1 = -2r \cos(2 \pi f_0 T) \\
      a_2 = r^2 \\
      f_0 = \text{resonant frequency} \\
      T = \text{sampling period} $

or more "refined:"

$$y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2]$$
where $ b_0 = \frac{1 - r^2}{2} \\
    b_1 = 0 \\
    b_2 = -b_0 $

The text says of the parameter $r$ that "the closer $r$ is to 1.0, the narrower the bandwidth of the resonance peak." This sounds to me like the definition of $Q$, but I suspect that it's not identical.

My question: What is the relationship between $r$ above and the more usual filter parameter $Q$?

I am just beginning to learn (self-teach) how to build filters, and I am still unclear on some of the concepts. Thanks in advance for any help.

See also this question and its answer.

I'll give you one quick answer (without proofs).

If your mapping from $s$ to $z$ is the Bilinear Transform, then

$$ Q = \frac{\sin(2 \pi f_0 T)}{2} \frac{1 + r^2}{1 - r^2} $$

If it's Impulse Invariant, I think it's

$$ Q = -\frac{\pi f_0 T}{\ln(r)} = -\frac{2 \pi f_0 T}{\ln(r^2)} $$

For high $Q$ and reasonably low $f_0 T$, the two expressions come out nearly equal.


Additional prefix and suffix required when windowing an OFDM symbol? To reduce out of band energy, it is common to apply a Raised Cosine filter to the front and back of an OFDM symbol. It is also common to overlap the windowed portions of adjacent symbols to eliminate amplitude dips.

For an OFDM symbol consisting of a Cyclic Prefix followed by NFFT points from an IFFT, is it common to insert an additional prefix and suffix for the windowed portion? If that isn't done then the CP and the IFFT data are shaped by the window, distorting the constellation. A pretty good description of what I am talking about is here:

http://zone.ni.com/reference/en-XX/help/373725A-01/wlangen/windowing/

Yes, it is necessary to add an additional pre/suffix to the OFDM symbol that corresponds to the window length when windowing is applied. As the length of the already existing cyclic prefix is usually chosen as the maximum delay spread of the channel no further Inter-symbol interference (ISI) must be introduced or otherwise orthogonality is lost. The windowing method described in the link you referenced is introducing additional ISI by overlapping consecutive symbols.

Furthermore, weighting some samples of the OFDM time-domain signal with different factors is a distortion that also distorts the constellation in frequency domain, as you say.

Note that windowing is not filtering. Windowing is the multiplication of the OFDM signal $x_k$ with some window function $w_k$. In contrast, a filter has a memory.


Digital implementation of first order analog filter using bilinear transformation I'm trying to create a digital filter from a first order analog filter with transfer function
$$H(s)=\frac{1}{1+\tau s}$$
with time constant $\tau=.1\text{s}$, and  sampling rate $f_s=1000\text{Hz}$.

Applying the bilinear transform in Matlab however appears to yield a filter with the a different 3dB point than expected. I expect the 3dB point to be at $\frac{1}{\tau}=10\text{Hz}$, but it appears to be around $1.6\text{Hz}$. Any idea what I could be doing incorrectly?

Matlab code:

fs = 1000;
tau = .1;

num = 1;
den = [tau, 1];
[numd,dend]=bilinear(num,den,fs);

[h, f] = freqz(numd,dend,4096, fs);

figure(1); clf();
subplot(211); semilogx(f,20*log10(abs(h))); hold on
plot([.1, 1000], [-3 -3],'r');
grid on; ylim([-40,1]); ylabel('gain (db)'); xlim([.1, fs/2]);
subplot(212); semilogx(f, angle(h)*180/pi);
grid on; ylabel('phase(rad)'); xlim([.1, fs/2]); xlabel('frequency(Hz)');

The expected 3dB frequency is wrong because of radian conversion. With $s=j\omega$,
$$H(j\omega) = \frac{1}{1+\tau j\omega} $$
and $20\log_{10}|H(j\omega)|\approx-3$
when $\omega=\omega_{3dB}=\frac{1}{\tau}$. However $\omega=2\pi f$, so
$$f_{3dB}=\frac{\omega_{3dB}}{2\pi}=\frac{1}{2\pi\tau}.$$
With $\tau=0.1\text{s}$, $\ f_{3dB}=1.59\text{Hz}$ as on the original plots.


How to get the conjugate symmetry of the complex matrix using MATLAB? How to get the conjugate symmetric of the complex matrix using matlab

This question is off-topic, but:

At = A';

At is the conjugate transpose of A.


How to design a lowpass filter using Hamming window and sinc function I have these three equations and have a noisy voice signal. I have to design a low pass filter of length N=1000 using hamming window technique. I have selected its cutoff frequecny 0.25 (normalized).
\begin{eqnarray}
h_{id}[n] &=& \frac{\sin{(2\pi f_c (n - \frac{N-1}{2}))}}{\pi (n - \frac{N-1}{2})},\ 0\le n  \le N-1 \\
  w[n] &=& 0.54-0.46\cos{(\frac{2\pi n}{N})} ,\ 0\le n \le N \\
 h_F[n] &=& h_{id}[n]\cdot w[n] \\
\end{eqnarray}
Now the problems that I am facing are: if I follow the equations, the dimensions of both $h_{id}[n]$ and $w[n]$ becomes unequal and I cannot multiply them to get the finite length impulse responce of sinc function; also as I am a newbie I really do not know how to apply this filter to the voice signal.

Note: According to the instructions from the Professor, I can not use any built-in function or commands except the "filter" command

Any help would be appreciated.

The Matlab code:

clc; clear all;  close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 1000;
n = 0:N-1;
new_n = 0:N;
fc = 0.3;
temp_N = (N-1)/2;

h_id = sin(2*pi*fc *(n - temp_N))./(pi *(n - temp_N));
figure('Name','Voice Signal in Time Domain','NumberTitle','off');
subplot(2,1,1)
plot(n,h_id)
title('Continous Time Plot (not true: as we have a DT signal, it is just for understandings)')
xlabel('Time \rightarrow')
ylabel('Magnitude \rightarrow')
axis([450 550 -0.2 0.7])
subplot(2,1,2)
stem(n,h_id)
title('Discrete Time Plot')
xlabel('Samples \rightarrow')
ylabel('Magnitude \rightarrow')
axis([450 550 -0.2 0.7])

w = 0.54 - 0.46*cos((2*pi*new_n)/N);
h_F = h_id .* w;
freqz(h_F)

fileName = 'D:\Studies\Semester#6\Digital Signal Processing\Assignments\Assignment 2\VOICE.wav';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[y, fs] = audioread(fileName);
N_FFT = length(y);
time = (1/fs) * N_FFT;
t = linspace(0, time, N_FFT);
figure('Name','Voice Signal in Time Domain','NumberTitle','off');
plot(t,y);
title('Time Domain plot of the Voice Signal')
xlabel('Time \rightarrow')
ylabel('Magnitude \rightarrow')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

playerObj = audioplayer(y, fs);
play(playerObj);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signalFFT = fft(y);
figure('Name','Voice Signal in Frequency Domain','NumberTitle','off');
subplot(2,1,1)
xaxis = (0:N_FFT-1)/N_FFT;
plot(xaxis, abs(signalFFT));
title('Frequency Domain Plot of the Voice signal (without shifting[fftshift])')
xlabel('Normalized Frequency \rightarrow')
ylabel('Magnitude \rightarrow')
subplot(2,1,2)
plot(xaxis, fftshift(abs(signalFFT)));
title('Frequency Domain Plot of the Voice signal (with shifting[fftshift])')
xlabel('Normalized Frequency \rightarrow')
ylabel('Magnitude \rightarrow')

If you can’t make the lengths of your Sinc kernel and window function the same, you can zero-pad the shorter of the two before you multiply them together.


Phase response of FFT in practice As I know if the impulse response is symmetric around sample zero phase response should be entirely zero.

The code below just set a rectangular window for vector "in"

const int N = 10;

// in = 0 , 0 , 0 , 0 , 1 , 1 , 1 , 0 , 0 , 0 , 0
std::vector< std::complex<double> > in (N);
std::vector< std::complex<double> > out (N);
std::vector< std::complex<double> > polarOut (N);

auto middleElem = in.begin() + in.size()/2;
std::fill( middleElem - 1, middleElem + 2, 1);

fftw_plan my_plan  = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(&in[0]),
                                         reinterpret_cast<fftw_complex*>(&out[0]), FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(my_plan);

std::transform( out.begin(), out.end(), polarOut.begin(),
                []( auto& in ){
                                return  std::complex<double>( std::abs(in), std::arg(in) );
                              } );

I getting back a sinc like result in magnitude,
But I am so curious why I have values different than zero in phase :

    polarOut    <11 items>  std::vector<std::complex<double>>
        [0] (3.000000, 0.000000)    std::complex<double>
        [1] (2.682507, -2.855993)   std::complex<double>
        [2] (1.830830, 0.571199)    std::complex<double>
        [3] (0.715370, -2.284795)   std::complex<double>
        [4] (0.309721, -1.999195)   std::complex<double>
        [5] (0.918986, 1.427997)    std::complex<double>
        [6] (0.918986, -1.427997)   std::complex<double>
        [7] (0.309721, 1.999195)    std::complex<double>
        [8] (0.715370, 2.284795)    std::complex<double>
        [9] (1.830830, -0.571199)   std::complex<double>
        [10]    (2.682507, 2.855993)    std::complex<double>

Any explanation about the phase response of behaviour of this fft will be great.

The problem is that your input vector is not symmetric when considering the definition of the DFT:

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

Note that from (1) the time indices start at $n=0$. If you periodically continue the signal $x[n]$ in your example you get

n:    ... -5 -4 -3 -2 -1 | 0 |  1  2  3  4  5  6  7  8  9  10

x[n]: ...  1  0  0  0  0 | 0 |  0  0  0  1  1  1  0  0  0  0

If $x[n]$ were symmetric, the values for negative $n$ would need to be the same as the values for positive $n$. For $n=\pm 4$ this is not the case. If you defined your signal instead as

x[n] = 0 0 0 0 1 1 1 0 0 0

i.e. a signal of length $N=10$, then it would be symmetric and, consequently, its DFT (FFT) would be real-valued, apart from numerical errors.


What Is the Relation Between Deblurring and Deconvolution in Computer Vision and Image Processing? The deblurring problem can be modelled as follows

$$
f = \phi u + \epsilon, \; \epsilon \sim N(0, \sigma)
$$

where $\phi$ is a filter (e.g. a low-pass filter) and $\epsilon$ is a Gaussian noise.

In computer vision, what is the relation between deblurring and deconvolution?

Let me present the following Diagram:

So, both Deblurring and Deconvolution are operations within the family of Image Restoration (Which is a subset of Inverse Problem set).

Basically we build the Image Restoration set by different Degradation Models.
The one related to the question are:

Linear Degradation Model
Namely, the degradation is made by a Linear Operator.
Spatially Invariant Model
A Model where the degradation is the same for any place in the image.
We created it as a subset of Linear Model though it doesn't have to be. But for clarity.
Any operator which is both Linear and Spatially Invariant can be defined by a Convolution Operation. Hence it can be reversed by Deconvolution.
Low Pass Operator
A set of images degraded by a Low Pass Filter. Namely by a convolution with a Low Pass Filter.

Now, there is a set of degradation which basically create a blurry image.
Reversing this operation is called Deblurring.
In case the blurring is made by a Low Pass Filter applied by Convolution the Deconvolution of this operation is also a Deblurring process.


Help with finishing this integral, to obtain the power spectral density of a pure cosine wave I am trying to evaluate the power spectral density $S_{xx}(f)$ of a cosine signal $x(t) = A\cos(2\pi f_0t)$, by starting from its definition for deterministic power signals
$$S_{xx}(f) = \lim_{T\rightarrow\infty}\frac{1}{T}\left|X_{T}(f) \right|^2\;, \tag{1}$$
where $X_T(f)$ is the Fourier transform of a signal $x(t)$ which has been truncated to a time window $T$, given by $$X_T(f) = \int_{-T/2}^{T/2}x(t)\;e^{-2\pi if t}\;dt \;.\tag{2}$$

I know that the answer should be $$S_{xx}(f) = \frac{A^2}{4}\bigg[\delta(f-f_o) + \delta(f+f_o)\bigg]\;,\tag{3}$$ but I am stuck with my working. Can anyone check it so far, and help me fill in the gaps? Perhaps my mathematics is not as good as it ought to be!

\begin{align}
  S_{xx}(f) &= \lim_{T\rightarrow\infty}\frac{1}{T}\left|X_{T}(f) \right|^2\ \tag{4}  \\
            &= \lim_{T\rightarrow\infty}\frac{1}{T}\left[ \int_{-T/2}^{T/2}x(t)\;e^{-2\pi if t}\;dt  \int_{-T/2}^{T/2}x^\ast(t')\;e^{2\pi if t'}\;dt' \right] \tag{5} \\
            &= \lim_{T\rightarrow\infty}\frac{1}{T}\left[ \int_{-T/2}^{T/2}\int_{-T/2}^{T/2} A^2\cos(2\pi f_0t)\cos(2\pi f_0t')\;e^{2\pi if(t'-t)}\;dt \;dt'\right] \tag{6}  \\
            &= \lim_{T\rightarrow\infty}\frac{A^2}{2T}\left[ \int_{-T/2}^{T/2}\int_{-T/2}^{T/2} \Big[\cos\big(2\pi f_0(t+t')\big) + \cos\big(2\pi f_0(t'-t)\big)\Big]\;e^{2\pi if(t'-t)}\;dt \;dt'\right] \tag{7}
\end{align}
where the last line is obtained by using the product-to-sum identity for cosines. Now, I'm a little stuck. I can make the substitution $\tau = t'-t$, which simplifes the exponential

$$  S_{xx}(f) = \lim_{T\rightarrow\infty}\frac{A^2}{2T}\left[ \int_{-T/2-t}^{T/2-t}\int_{-T/2}^{T/2} \Big[\cos\big(2\pi f_0(2t+\tau)\big) + \cos\big(2\pi f_0\tau\big)\Big]\;e^{2\pi if\tau}\;dt \;d\tau\right] \tag{8}  $$

but now I am not sure how to proceed. I can also replace my cosines with complex exponentials to obtain

$$\lim_{T\rightarrow\infty}\frac{A^2}{4T}\left[ \int_{-T/2}^{T/2}\int_{-T/2}^{T/2} \Big[ e^{2\pi if_0(t+t')} + e^{2\pi if_0(t-t')} +e^{2\pi if_0(t'-t)}+e^{-2\pi if_0(t+t')}    \Big]\;e^{2\pi if(t'-t)}\;dt \;dt'\right]$$

Can someone teach me how to finish this? Ideally I would like it if someone could show me how to evaluate it in terms of $T$, and then show how the $\delta$ functions arise as we finally take the limit $T\rightarrow\infty$.

You don't need to make it so complicated as to use two variables $t$ and $t^\prime$.

\begin{align} X_T(f) &= \int_{-T/2}^{T/2} x(t)\exp(-j2\pi ft)\,\mathrm dt\\
&= \int_{-T/2}^{T/2}\frac{\exp(j2\pi f_0t)+\exp(-j2\pi f_0t)}{2}\exp(-j2\pi ft)\,\mathrm dt\\
&= \frac 12\int_{-T/2}^{T/2} \exp(j2\pi (f_0-f)t)+\exp(-j2\pi (f_0+f)t)\,\mathrm dt\\
&= \frac 12 \left[\frac{\exp(j2\pi (f_0-f)t)}{j2\pi(f_0-f)}~ + ~\frac{\exp(j2\pi (f_0+f)t)}{j2\pi(f_0+f)}\right\vert_{-T/2}^{T/2}\\
&= \frac 12 \left[\frac{\exp(j\pi (f_0-f)T)-\exp(-j\pi (f_0-f)T)}{j2\pi(f_0-f)}\right.\\
&\qquad\qquad + \left.\frac{\exp(j\pi (f_0+f)T)-\exp(-j\pi (f_0+f)T)}{j2\pi(f_0-f)}\right]\\
&= \frac T2 \big(\operatorname{sinc}((f_0+f)T) ~+ ~\operatorname{sinc}((f_0-f)T)\big)\\
&= \frac T2 \big(\operatorname{sinc}((f-f_0)T) ~+ ~\operatorname{sinc}((f+f_0)T)\big).
\end{align}
So we get \begin{align}\frac 1T\big|X_T(f)\big|^2 &= \frac T4 \big(\operatorname{sinc}^2((f-f_0)T) ~+ ~(\operatorname{sinc}^2((f+f_0)T)\\
&\qquad\qquad + 2\operatorname{sinc}((f-f_0)T)(\operatorname{sinc}((f+f_0)T)\big)
\end{align}
In the frequency domain, the two $\operatorname{sinc}^2$ functions are of height $\frac T4$, are centered at $\pm f_0$ with central lobes of width $\frac 2T$ Hz. As $T \to \infty$, the central lobe widths shrink to $0$ and it is straightforward to show that for any value of $f$ other than $\pm f_0$,
$$\lim_{T\to\infty} \frac 1T\big|X_T(f)\big|^2 =0, ~ f \neq \pm f_0.$$
Now, $\frac 1T|X_T(\pm f_0)|^2$ diverges to $\infty$ and with the usual trick of pulling a Dirac delta out of the hat in such cases, we get that it must be that
$$S_{xx}(f) = \lim_{T\to\infty} \frac 1T|X_T(f)|^2 = \frac 14 \big(\delta(f-f_0) + \delta(f+f_0)\big), ~~ \text{when} ~x(t) = \cos(2\pi f_0 t).$$


How to get the data (pointwise frequency, amplitude and phase) in a FFT performed on ImageJ? How to find pointwise readouts of the amplitude, frequency and phase of the DFT underpinning the FFT image?

Once an image is FFT-ed in ImageJ, placing the cursor over any points on the FFT plot results in an output of values at the bottom of the main menu, as in this example:

According to the documentation it simply outputs the location of the cursor:

  If the mouse is over an active frequency domain (FFT) window, its
  location is displayed in polar coordinates. The angle is expressed in
  degrees, while the radius is expressed in pixels per cycle (p/c). The
  radius is expressed in [units] per cycle (e.g. mm/c) if the spatial
  scale of the image was defined using Analyze/Set Scale.

I don't know if there is any additional information within this line to link to the actual sine (or cosine) wave at this particular point in the Fourier plot:

r = .18 p/c (644), theta = 55.12 deg, value = 121

Some playing with ImageJ has given me some preliminary intuition of what these numbers signify. Picking up extreme positions of the cursor as it hovers over the FFT image or plot (yellow cross on the following images) makes it easier to interpret. First, placing the cursor very close to the bright dot in the center and along the x axis yields the following read-out:

This is the slowest frequency in the image, and corresponds to 1 cycle from side-to-side. The image is 256 x 256, and hence, the frequency is 256 pixels / cycle.

What happens if we move the cursor a tiny bit further away to the left along the x axis:

Right, now the frequency has doubled to 4 cycles from side-to-side, which corresponds to 64 pixels/cycle, which (given enough patience and skepticism) can be proven by magnifying and counting the pixels:

What happens at the Nyquist limit:

The limit is 2 pixels/cycle for 128 cycles.

The direction, naturally can be at any angle around the clock face:

and if we select a darker pixel (as in this case), the value will consequently be low.

It gets interesting at the corners, where the PT tells us that the maximum number of cycles will be $\frac{\sqrt{256^2+256^2}}{2}=181$:

There ought to be an easier method, but awaiting a formal answer this would possibly get the job done:

The read-out of the information captured on the image may go like this:

At the point with the cross hairs on the frequency space:

  $\text{Freq.} = 2.245$ pixels/cycle - on a $256^2$ image, $256/2.245 = 114$ cycles in the direction of...

  $\text{Orientation} = \theta = 48.2^\circ$

  $\text{Amplitude} = \vert X[k]\vert=\sqrt{X_\text{Re}^2+X_\text{Im}^2}=\sqrt{-244.868^2+456.238^2}=384.95$

  $\text{Phase} =\arctan\left(\frac{X_\text{Im}}{X_\text{Re}}\right)=\arctan\left(\frac{456.238}{-244.868}\right)=-61.78^\circ$

whether this is correct, and whether there is an easier way to get this info without having to calculate on the side these values, or having to use several separate, and disjointed (FFT and Complex) output windows, justifies keeping the post open.

Here is a very practical illustration. The credit goes to @Herbie and @emartini.

NOTE: The complex window has a slide bar at the bottom, which has to be drawn to the left. Also, it is necessary to click twice on the FFT window without moving the cursor to get the three output values.

The output is concordant with Matlab with adjustments nicely described by @SleuthEye here and here.

Here is an example:

with ImageJ:

Frequency = 10.24 pixels/cycle (25 cycles)
Theta (direction) = 16.26 degrees
Real part = -1.255
Imaginary part = 10.142
Phase = arctan(10.142 / -1.255) = -82.95 degrees
Magnitude = sqrt(10.142^2 + 1.255^2) = 10.2194

and with Matlab:

Im = imread('circ.png');
pkg load image
Im = rgb2gray(Im);

A = fft2(double(Im) / 255);
Ashifted = fftshift(A);
Ashifted(121,153)

i = 121;
j = 153;

center = size(A) / 2 + 1;
dx = (j - center(2)) / size(A,2);
dy = (center(1) - i - 1) / size(A,1);

direction = (atan2(dy, dx))
dir_degrees = direction * (360 / (2*pi))
frequency = 1 /sqrt(dx*dx + dy*dy)

The output:

ans =  -1.2553 + 10.1425i
direction =  0.28379
dir_degrees =  16.260
frequency =  10.240


(Newbie question) Calculating the frequency spectrum I'm quite new to this topic, with some basic knowledge however.
I need to analyse an audio signal for different characteristics.
Like the audio editor I use (an older Version of Audition, v.3)
I would like to calculate the occuring frequencies.

For every time stamp and
the overall occurences of the frequencies of the whole audio.

If I'm not totally wrong then for calculating the frequency spectrum
I have to take an interval (window) of samples.
The DCT/ FFT of that returns the complex fourier coefficients
(whose magnitudes are the amplitudes of the frequencies).
Correct so far?
What I do not understand:
What is the frequency spectrum of a specific time stamp then?
Or does it only make sense for a complete interval/ window of samples?
But Audition shows a different frequency spectrum for every time stamp...
Is it something like a "sliding window" there, calculated for the center sample?
Associated with the above:
If it should be (the middle of a) sliding window:
How to calculate spectrums for near-border timestamps
(which don't have complete surroundings)?
And how to calculate the overall occurences of the frequencies
of the whole audido?
I can't use the sliding windows for that, as it would count
frequencies multiple times (due to the overlapping) ...
Is it the average of all the side-by-side (= non-sliding) windows?
Or am I totally wrong with it?
Again: I hope my questions don't hurt too much. ;-)

I'm almost tempted to answer with a counter question: what is frequency? It's a question I like to ask students that have already taken signals and systems courses and it's surprising how often they don't know what to answer.
While the whole answer involves some more discussion, what we can definitely say is that frequency is associated to a rate of change of the signal. That automatically means that looking at a signal in one moment of time makes it impossible to say anything about its frequency content, because we cannot see how rapidly it changes.
Therefore, most of what you say is true: a frequency analysis only makes sense for a certain window in time. One can use sliding windows, typically one defines some overlap between adjacent windows, depending on what exactly one is after. The length of the window is a trade-off: longer windows will give you a better frequency resolution (watching the signal change over a longer period of time makes us more certain how exactly it varies) but less resolution in time (if the frequency content changes over time, we may want to be sensitive to that change and if a change is so sudden it appears within one window, we might miss it).
For the borders, it depends what you know about your signal. Either you clip them, or you pad them with something that is meaningful. This can be zeros, it can be a mirrored image of the samples you know (this is more common in image processing), it can be a periodic copy of the end (this is done in communications for other reasons). In most cases, signals are long enough so that we don't need to worry about the edges that much.
As for your last question, I'm not sure I get what you mean by overall occurences. But of course you could average your windowed spectra to get an "average spectrum". This is referred to as Bartlett's method, some sources call this the periodogram (careful, the term is used in conflicting manners in the body of literature).


Do I have to flip my kernel when performing an FFT-based convolution? I have a 2D image and I convolve it with a 2D kernel image using FFT. So far I was always using symmetric kernels (e.g., Gaussian with stddev_x = stddev_y). As a result, I never bothered thinking whether I have to flip my kernel image or not because it wouldn't have made any difference. However, now I want to convolve my image using an elliptical Gaussian kernel with stddev_x != stddev_y and an arbitrary angle. Now, do I have to flip my kernel image prior the FFT convolution? Or the flipping is required only when using the usual convolution algorithm and not the FFT-based one?

Thank you.

Flipping & dragging is an animative method used in graphical computation of convolution in time (or space) domain. It is the result of an argument manipulation in $h[n-k]$ (or $x[n-k]$) signal animated as a function of n, but drawn on an axis of k, in the convolution sum: $$y[n] = x[n]*h[n]=\sum {x[k]h[n-k] }$$ On the other hand in frequency domain DFT (FFT) based convolution implementations you don't need it as you simply multiply the two DFTs.


Nonlinear tricks for blocking wide band noise Conventional methods of filtering allow one to block parts of a frequency band while passing others.

But thinking outside of conventional means, and perhaps not based on any frequency methods at all, are there filtering methods for blocking wide band noise, however passing rapid moving signal?

For example. Say I have a square wave with 'hair' on it (white-like noise). I just want to get the square wave with corners as sharply as possible - and minimum hair on the output. Outside of conventional FIR, IIR, Kalman filtering are there any nonlinear tricks?

A very common non-linear filtering technique for removing noise while preserving edges is median filtering. It basically replaces the samples inside a window by their median. Total variation denoising is another non-linear technique that works especially well for denoising piece-wise constant signals. You can find some more information about this method and Matlab software here.


Estimating Low Frequencies (including amplitude and phase) For the past month I've been chasing the DFT/FFT rabbit.
Along the way I've been cut-n-pasting the name of other rabbits to chase later:
Autocorrelation, Phase Vocoder, Parabolic Approximation, Wavelet Transform

Here's my problem.
When using an FFT, the lower frequencies are always going to be missing.
I actually found the same problem described here:
Problems with the FFT/IFFT

Frankly it is worse than that. The energy of these frequencies are not missing at all.
They are dissolved into the first several bins, making these bins semi-useless for interpolation.

For now, I simply want a general approximation of the Low Frequencies, including amplitude and phase. In FFT terms, I want bin #0.05, 0.10, 0.15, 0.20, 0.25 ...up to say around bin 3.50. Realistically, I understand that anything to small, say 0.05, would almost amount to Linear Regression. But in my mind, 0.5 or 0.75 should be doable.

I even toyed with changing k to a double (which was 'cough' "interesting"):
best window for close frequency components

When Moses came down from the mountain, he did not carry tablets of stone on which it was commanded that in the equation: F = k*(df) = k*sample_rate/N, that k has to be an integer.

Anyway, I've decided it is time to try something else.
So, what rabbit should I chase next?
(even if it is close approximation)

PS: recommending a bigger sample is not the answer - I'm already sampling my available limit.

Edited to add pics (hoping for a possible new answer):

Test Case (I want a Blue approximation, even if lousy):

Rectangle (of course, lousy Green and so-so Red):

Gaussian 3.0 (only Red can be Interpolated, but decent):

Smaller Gaussian gives a lousy Green Interpolation and decent Red.

Edited to answer general questions (I'm sorry I did not do this to begin with):

What's the mathematical description of the data in the plot?
I'm guessing Sum of Sinusoids (with noise).

What's the sample rate?
All of my samples are at different rates, as in, each entire file size is used.
Rather than Hertz/periods/cycles, the measurements are in deserialized data units.
Most often I work with 8-byte doubles, but sometimes 4-byte integers.
Therefore, with doubles, you can say that my Fs = fileBytes/8

How many samples can you capture?
One each per file.

What's the application?
I have no clue - it's for a client.

In the frequency representation given by the FFT there are no frequencies missing and their is no binning of frequencies.

FFT decomposes a discrete signal in to a sum of sinusoids of certain frequency and phase. This representation is exact and has the same energy as the original signal. No information is lost, as you can do an inverse FFT and recover the original signal exactly. The binning is coming from the discretisation of the signal, not the FFT.

Also for a finite length signal, there are no frequencies lower than the one with the period equal to the length of the signal. This is from the mathematical definition of the spectrum. Put more simply, the DFT/FFT 'assumes' the signal is repeating. Zero padding of the signal therefore increases the period and therefore the lowers the first frequency.

If there are lower frequencies due to some prior knowledge, e.g. if we sampled a sine wave over half its period, then this information must be included.

One solution may be to roll your own spectrum and scale the frequency components with period longer than that of the signal. The assumption being that the low frequency component continues.

For example, let $f_S$ be the frequency you are interested in, with period = $1/f_s > 1$ longer than the period of the signal which is normalised to 1.

Now you convolve the signal with the sine and cosine waves at that frequency. No zero padding.

$$f_e = \sum_t^N f[t] \times cos[2\pi f_s t]$$
$$f_o = \sum_t^N f[t] \times sin[2\pi f_s t]$$

Then multiply by $1/ f_s$ to account for the assumption the signal outside of the sample continues like this, and multiply by 1/N to account for number of samples.

$$f_e' = f_e \times 1/ f_s \times 1/N$$
$$f_o' = f_o \times 1/ f_s \times 1/N$$

Amplitude and phase are then given in the usual way.

This can be done at whatever frequency graduations you desire.

Edit:
Changes with regard to first comment.

Edit 2:
Answer riddled with errors. Updated.

Edit 3: Add factoring for N

--
Edits welcome if I have this wrong.


How to detect S wave peaks in ecg While there are plenty of methods to detect R peaks it appears that detecting S peaks is less discussed. Is there a recommended way of detecting S peaks?

Additionally, could you recommend an implementation in python?

If you are using Python you could try this:

def S_point(signal, R_peaks):

    num_peak=R_peaks.shape[0]
    S_point_list=list()
    for index in range(num_peak):
        i=R_peaks[index]
        cnt=i
        if cnt+1>=signal.shape[0]:
            break
        while signal[cnt]>signal[cnt+1]:
            cnt+=1
            if cnt>=signal.shape[0]:
                break
        S_point_list.append(cnt)
    return np.asarray(S_point_list)

If you allready have the R-peaks then you can use the signal and the R-peaks to get the S_points with this function

S_points=S_point(signal,R_peaks)

And finaly plot it like this

plt.plot(signal)
plt.plot(R_peak,signal[R_peak], 'yo')
plt.plot(S_points,signal[S_points], 'ro')


Techniques to reject noisy neural network input Suppose an artificial neural network is used to approximate a sine wave (shown in red in the graph below), given the linear input variable $x$ (scaled such that the ANN input is $x_{\rm nn}\in[-1;1]$).

The mean squared error between the target function (sine wave) and the trained ANN output is as shown approximately $4.8\times10^{-3}$, if the ANN structure is 1-8-1 (one input $x_{\rm nn}$, 8 hidden nodes and 1 output node).

Now suppose the ANN structure is changed to 2-8-1, where the second input is normally distributed random noise. The output is then as shown below. The MSE has increased to approximately $1.2\times10^{-2}$ (using the exact same noise sequence as during training) or slightly larger if a new noise sequence is used when testing the network.

The results above were obtained using standard gradient-descent type training with regularisation. All nodes in the hidden layer have $\rm tanh()$ activation functions and the output node is linear, as is standard with function approximation.

What methods are available to automatically remove the second input?
That is, training methods allowing the ANN to automatically discover the fact that the second input is noise, and therefore not a useful input. (Therefore excluding all methods where the human programmer must train several different ANN architectures with different input combinations, and evaluate all of them to decide which inputs to use and which to exclude.)

This question extends to larger ANNs where it is difficult to see beforehand which inputs are useful. Ideally, all weights associated with nonsense inputs should automatically be forced to zero.

It seems like your model is overfitting because you are not providing it with enough training data.

If I understood correctly, you want the neural network to learn $\hat{y} = \hat{f}(x,n)$ where $n \sim \mathcal{N}(\mu,\sigma)$ is noise and $x$ is in the interval $[-1,1]$ by minimizing
$$
\frac{1}{N}\sum_{i=1}^N |y_i - \hat{y_i}|^2 = \frac{1}{N}\sum_{i=1}^N |f(x_i,n_i) - \hat{f}(x_i,n_i)|^2
$$

And the function you are trying to learn is $f(x,n) = sin(\pi x)$, right?.

The network should learn to ignore the noise $n$ given that you provide it enough examples for training. That is, for different values of $x$, your training data should have a lot of examples with different values for $n$, all having the same target $sin(\pi x)$. If you show the network enough examples where the $n$ value is not useful for making the prediction, it should learn to minimize its effect on the output. The more training data you have, the better the network can learn to do this.


Sound quality with really high sampling rate, but really low bit depth with dithering? This question is probably more for fun and theory, but how will the sound quality be, if you create a sound with really high sampling rate (like 192kHz) and really low bit depth (e.g. 3bit)?

Can you emulate 44.1kHz/16bit with 192kHz/3bit (+dithering)?

I can probably think this through and post an answer myself with, sound examples computed in MATLAB, but maybe someone has already done that :)

well, there is a sorta information-theory theorem regarding this (called the Gerzon-Craven limit) that was published long ago by Michael Gerzon and Peter Craven back in the late 80s.  i can't get a direct link to the paper reference, but type in "Optimal Noise Shaping and Dither of Digital Signals" into the AES Convention paper site.  this paper is a quarter-century old.  i just discovered that i have a camera copy of that old AES preprint and can send it to you if you send me an email to rbj@audioimagination.com .

i also have a copy of a paper making reference to the Gerzon-Craven limit On psychoacoustic noise shaping for audio requantization .

there is an interesting title regarding this called Approaching the Gerzon-Craven Noise Shaping Limit Using Semi-Infinite Programming Techniques. i don't have a copy of that and do not know how to get a free copy.

essentially, because of the Shannon Channel Capacity theorem, the answer to your question is "yes almost, theoretically", as long as you don't care about the magnitude of the quantization noise outside of your original 22.05 kHz baseband.  this is because the bit rate at 16-bit, 44.1 kHz is about the same (maybe just a teeny bit more) as the bit rate at 3-bit, 192 kHz and you can make an information-theory case that, because the information rate of the two cases are about the same, then the S/N ratio in the band of interest is also about the same.

but the real answer is "no, practically", because of the difficulty in constructing a 3-bit noise-shaped modulator with such extreme specs and that remains stable.  with a higher sampling rate than 192 kHz, the specs get less extreme.

if you have a working system (what's it coded in, MATLAB?), i would be interested in your noise-shaped modulator.  what order it is and how well it behaves.  sometimes these nasty high-order ΣΔ modulators behave badly in the presence of silence and we get, what are called, "idle tones".  if you think you have a working solution, i would be interested in hearing how it works for very low-level input.


How are video thumbnails chosen? Of the thousands to millions of frames that a video might have, how is a particular representative frame chosen to be the thumbnail?

There are many ways to skin this cat. I'll explain one way that I know is used in industry.

Create candidate frames algorithmically using the image features; e.g., the histogram, sharpness, presence of humans, faces, and anything else you can think of.
Serve different candidates to different users and see which ones are popular
You can stop once the popularity contest results stabilize. However, for bonus points you can learn the ranking function that yielded the list based on the features enumerated in part one.
Use the ranking function to select the best frame, or cull the list of candidates for subsequent tests.


Adaptive filter with two inputs I have general theoretical questions:

Is it true that an adaptive filter with two inputs (one normal and one delayed by the single time increment) can completely get rid of any single frequency noise?
Is it then true that a three-input adaptive filter can (completely) get rid of a noise consisting of two harmonics, etc.?

Maybe there is a theorem for this, any pointers would be appreciated.

I've had a look at the problem in the book you've cited (Figure 10.7, Adaptive Filter for Noise Cancellation).  As Dan has stated, a linear combination can filter out a specific frequency, but it depends on the specific noise frequency, the sample rate, and the length of the adaptive filter delay-line (i.e. the number taps or weights used in the adaptive filter).

In general, however, a given adaptive filter arrangement used for noise cancellation (i.e. sample rate, length of delay-line, adaptive algorithm) will dictate what specific frequencies of the noise will be attenuated.


Why does the high frequency component of an image contain more noise? I am working with wavelet transforms. When we apply a discrete wavelet transform to an image, it decomposes into four components named the approximation component, thevertical component, the horizontal component, and the diagonal component.

These components are indicated as LL, HL, LH, HH respectively. I understand how an image is decomposed by using a wavelet transform, but I don't know why generally in an image high frequency component contain more noise.

There are different types of noise sources in imaging process.

Photoelectronic

photon noise: Poisson distribution, whole spectra;

thermal noise: Gaussian distribution, white noise, whole spectra;

Impulse

salt and pepper noise: mainly high frequency noise; (@user7358 made a very good point on the reason).

line drop: part or all of a line in the image lost, can be viewed as a deterministic high frequency noise;

Structured

periodic: low frequency noise;

aperiodic: JPEG noise, low frequency too;

detector striping: calibration differences among individual scanning detectors, mainly low frequency;

If you only observe the evident noise when you zoom in the high frequency band with wavelet, I assume it is mainly the shot noise (salt and pepper). When you are doing some simulations with noise added, try (1) salt and pepper,  (2) Gaussian, and (3) Gaussian with a high-pass band applied, and observe the difference after wavelet decomposition.


How poles are related to frequency response I have recently fallen into fallacy, considering pole s=1 as there is infinite response at frequency 1. Yet, response was only 1. Now, can you derive the frequency response, given the poles?

Secondly, the theory says that a system is stable when poles are in left s-plane and, thus, decay in time. But, wait. Does'n "pole" mean the infinite response -- the growth in time?

Finally, is it right question in DSP? IMO, D stands for digital whereas s-domain is analog. I do not find s-plane or Laplace transform tags to label my post.

update Thanks for the answers. It seems that I have got it except the one minor but fundamental thing -- the relationship of poles (and zeroes) with frequency. Basically, why are eigenvalues (or, how do you call the $s$ operator/variable) related with frequency? It should be somehow related with exponential growth and Laplace transform. I quite understand that poles happen to be eigenvalues (especially for discrete recurrences). But, how is this related with frequency?

One thing that really helped me understand poles and zeros is to visualize them as amplitude surfaces.  Several of these plots can be found in A Filter Primer.  Some notes:

It's probably easier to learn the analog S plane first, and after you understand it, then learn how the digital Z plane works.
A zero is a point at which the gain of the transfer function is zero.
A pole is a point at which the gain of the transfer function is infinite.
Often there are zeros or poles at infinity, which aren't always included in descriptions of the transfer function, but are necessary to understand it.
The frequency response in the S plane happens along the jω axis only.

The origin is 0 Hz, or DC, and the cutoff frequency of filters increases radially away from the origin.  Putting a pole at any point along a circle at a certain distance from the origin will produce the same cutoff frequency.
To increase the cutoff frequency of a filter, move the poles radially outward.
To increase the Q of a biquad filter, move the poles along the circle towards the jω axis, which keeps the cutoff frequency constant, but increases the effect that the pole has on the frequency response, making it more "peaky".

If a zero appears on the jω axis, then the frequency response will drop to zero at that frequency; if you input a sine wave at that frequency, the output will be 0.
If a pole appears on the jω axis, then the impulse response is an oscillator; any impulse will cause it to ring forever at that frequency.  Impulses have finite energy, but the response of the filter has infinite energy, so it has infinite gain.

A simple example is an integrator H(s) = 1/s:

This function equals 0 when s is infinite, so it has a zero at infinity.
This function equals infinity when s is zero, so it has a pole at zero.

In other words, it has infinite gain at DC (the step response of an integrator is forever-increasing), and the gain decreases as frequency increases:

Moving the pole away from the origin, along the imaginary axis into the left hand of the S plane, makes the gain at 0 Hz on the jw axis finite again, and now you have a low-pass filter:


About designing digital filters I'm currently using MATLAB's fdatool for filter design. Using that tool, I can easily design different kind of filters. For example, let's take a band-pass FIR filter with 10-40 Hz passband, and 5-10 Hz and 40-45 Hz transition bands. Usually, I design the filter with the selection "least-squares", which, if I understand correctly, uses the aforementioned method to find the best impulse response according to filter spesifications. To actually filter the signal, I use the command filtfilt, which does zero-phase FIR filtering.

Now, an alternative way to implement the filter would be to take the FFT of my signal, set frequencies outside the range 10-40 Hz as zeros, and then take the IDFT.

Is there any practical/theoretical difference between these two approaches? Will the frequency responses (magnitude and phase) be the same?

Filtering can be done in the frequency domain which is actually a very efficient technique (and it can very well be, that Matlab does this internally). However, for very long signals it's not as straight-forward as "taking the FFT and applying a frequency response". You can read up on overlap-add filtering which is a possibility to filter in frequency domain. However, this is mainly done because it can be faster.

However, this question is independent from filter design which is a completely different thing. Every FIR filter has a corresponding frequency response and it does not matter whether you convolve in time domain or apply the frequency response in frequency domain. So your approach becomes a question of filter design and the consequences it has.

Just setting unwanted frequencies to zero might completely eliminate those frequencies but it usually comes at the price of a significant ringing in the time domain which is usually unwanted (Wikipedia on ringing). So in fact, what you are proposing is just one way to design a filter and frankly not a good one in most cases.

By allowing a transition band in which the frequency response can gradually go from passband to stopband, degrees of freedom are gained that can be used to improve other properties of the resulting filter (for example eliminate the ringing or obtain a shorter impulse response). That's why Matlab implements so many different filter types, they all have different properties and selecting the most suited one is actually part of designing your signal processing system.

This topic is actually quite complicated and I suggest to read up on filter design. Whole books can be filled with this.


Filter response is changed according to its order I want to build a band-reject filter for the 60 Hz noise. So, I tried to use Butterworth filter to implement this specification. However, when I chose an order of 8, I got the following response, which makes sense.

However when I chose an order of 30. I got this weird response. Why is that?

The code to get the coefficient for the 8 and 30 order respectively is shown below

# fs = 200
b, a = signal.butter(8, [58 / (fs / 2), 62 / (fs / 2)], btype='bandstop')
b, a = signal.butter(30, [58 / (fs / 2), 62 / (fs / 2)], btype='bandstop')

An order of 30 is way too high, especially for a band reject filter with such a narrow band as you want.

Your filter will have 15 poles in a very small semicircle between frequencies 58 and 62 Hz (and another 15 poles in the negative frequencies).

Your a coefficients correspond to an order 30 polynomial with the 30 roots corresponding to the poles of the filter.

Note that any minuscule numerical approximation (which you are not free from) will alter the position of all roots. You will very probably get roots ending up outside of the unit circle (which makes your filter unstable).

So the ugly frequency response for the high order filter is expected.

As a rule of thumb, do not try to implement IIR filters of order above 10.


Phase information from product of complex conjugate transfer functions For an input signal $x(t)$ and output signal $y(t)$ through an LTI system I $H(t)$ I have found and interesting property that concerns signals' power spectral densities:

$${\lvert H(\omega)\rvert}^2 = \frac{\phi_y(\omega)}{\phi_x(\omega)}$$

This relation gives me information on the magnitude of $H(\omega)$ but no phase information. Is there a way to obtain phase information of $H(\omega)$ from this relation?

The answer to your question is no.

The relationship you show is well known and used to relate the input to output of Gaussian stochastic processes through a linear filter. Any filter that has the same magnitude response $\mid H(\omega ) \mid$ will produce the same result.

One can often infer a phase by assumptions on the filter such as linear phase or minimum phase but this isn’t  usually the case.

If you want the actual phase, you need to compute the cross spectrum.


How to model airborne sound channel I'm making a soft demodulation/decoding to communicate a microphone and a speaker in short range. I have to calculate the LLR using my constellation and the airborne sound channel probability density function (pdf). I have searched a lot about airborne sound or acoustic channel model, and I just found out some considerations about the attenuation depending on the frequency. My questions are:

Wich book can I read, so i can learn more about soft demodulation examples: BPSK, QPSK, QAM, FSK?
If I use a pdf estimator of the channel, how can I use the estimation to find the LLR of my received symbols?
Can I use another model to aproximate the airborne acoustic channel (it could be related to AWGN), or which book can I read to learn about how to model a communication channel?

Well, your channel model might give you an explicit PDF, but:

You need to realize that a channel is not a scalar, ie. your PDF doesn't describe a random variable $\in \mathbb C$ or $\in \mathbb R$, but usually something like an impulse response, that is, something like a function of delay (so, your random variable is $\in \mathbb R \times \mathbb R$ or $\in \mathbb C \times \mathbb R$).

I'm not an acoustics guy myself, but if I've learned one thing it's that at least for indoor acoustics, the "flat channel assumption" (that lets us wireless communication folks model channels as simple complex number) doesn't hold, and that your channel is very frequency-selective. Often, it's time-variant, and sometimes even non-linear.

So, I'm not quite sure what to recommend here:

If you're looking for soft decision in wireless communications, I personally think that Proakis is quite OK. That's one of the standard digital communications textbooks.
That's up to your channel model, and you'll have to specify that.
An AWGN channel is certainly not an adequate approximation of an audio channel.


LPC Formant Estimation: Why resample? The Burg algorithm implemented in PRAAT for formant estimation also works via LPC. First the sound is resampled to a sampling frequency of twice the value of maximum formant. Then a pre-emphasis is applied and finally, PRAAT applies a Gaussian like window to compute linear prediction coefficients through the Burg algorithm.

I understand that pre-emphasis is helpful for formant extraction by LPC since the human voice is not spectrally flat, but am confused about why Praat is resampling.
https://github.com/praat/praat/blob/bac77ad76481c65703ade8d7ddb190e9d67048d6/fon/Sound_to_Formant.cpp#L381

You want to resample down to alias frequencies higher than your maximum desired formant.
If you understand aliasing in the z-transform and the transfer function of a system then this should be obvious.