Does Increasing Resolution Lead to Have a Poor Resolution? As it has proven in some paper and books the time and frequency resolution of Morlet wavelet is :
  $\Delta t=\frac{f_c\sqrt{f_b}}{2f_i}$
 $\Delta f=\frac{1}{2\pi f_c \sqrt{f_b}}$
I don't know how these relations have obtained, I have asked it in detail in here.you can answer that if you want.
 but lets accept these relationship for now. there is some misunderstanding for me in here about the concept of resolution.
for this case, by increasing the term $f_c\sqrt{f_b}$ we have finer frequency resolution but poor time resolution. I am sure about finer frequency and poor time resolution because my practical results meet it.  but if you look at $\Delta t$ and  $\Delta f$ relationships you receive that increasing the term $f_c\sqrt{f_b}$ leads to decrease $\Delta f$  and increase the $\Delta t$.
does it mean that smaller value of resolution spell to have finer resolution and vice verca?

Indeed, that's the Heisenberg Uncertainty Principle - you can't have both very good frequency and time resolution. You always have to sacrifice something. In case of Short Time Fourier Transform it's straightforward, but for wavelets are being 'squeezed', which is changing their frequency resolution. Figure below describes more than a thousand words:

EDIT:
Because you are still trying to understand why it is like that, then please mind that by increasing the frequency of the wavelet ($f_c$), you should also shrink it in time (decrease $f_b$). For more info please refer for example to: Paul Addison - The Illustrated Wavelet Transform Handbook. You can find there a following graphs (Section 2.12):


Larger image after wavelet reconstruction Depending on the size of the original image and depending on the kind of padding used during 2D discrete wavelet transform and reconstruction, the reconstructed image has a larger size.

Is there a general method given a mxn matrix and given some padding mode to crop the obtained pxq matrix to the mxn matrix (without trying every possibility in order to find the one with the least error)?

A generic method consists in adapting the wavelet decomposition on the border, ie the wavelet varies on shorter frames. It can be called "wavelet on the interval" or be addressed with lifting schemes. But this requires caution and complicated computations sometimes.

A more lazy option (if I remember well)  is to use a pair of symmetric/anti-symmetric $M$-band wavelets, over $J$ levels. If you expand the image to a size that can be divided by $J^M$ (the most symmetrically as you can), and use  a symmetric extension that may reveal non-expansive if properly tuned. This may vary with implementations.


Sensing matrix for compressed sensing What are the differences between random binary sensing matrix  and
random Gaussian sensing matrix?
What the advantages and disadvantages of each matrix?

How can I choose the suitable matrix for a certain signal?

A sensing matrix maps input vector to measurement vector through linear wighted summation of input. What makes a specific matrix good, is application dependent. Now, both distributions more or less satisfy RIP. However hardware implementation of the Bernoulli matrix (binary or bipolar) is much much easier especially in analog domain. A Bernoulli wight is either 0 or 1 (or -1/1 in case of polar Bernoulli), but a Gaussian wight is a floating point figure. Multiplication of a flouting point number either in digital or analog, is resource consuming, while multiplication of a Bernoulli wight is feasible through implementation of a simple switch in analog domain or an and instruction in digital. As an example consider RMPI analog to information devices which compressively sample the signal in analog and then reconstruct the signal in digital. Prior to quantization the sensing matrix should be applied, were through incorporating a Bernoulli matrix, the multipliers are implemented as simples switches.

I also think, Gaussian Matrix has tighter RIP than Bernoulli (especially binary Bernoulli). As an example assume we have an sparse matrix  $X=[0,1,0,2,0]^T$ , and assume we are measuring it using a random matrix with two rows and 5 columns. A binary Bernoulli, is more probable to multiply one of the non sparse entries by zero and hence miss its contribution in the measurements vector.
Also, having only 1,0 and -1, the matrix rows are more probable to be similar and hence more coherency occurs.
Although it seems not very practical, in this paper the authors presented an analog circuit for application of Gaussian Matrix.


Amplitude value after a Discrete-Time integration calculation? I am using Matlab block Discrete-Time
Integrator on Simulink and I would like to know in advance, what  output of that block regarding to its amplitude and the current sampling frequency would be.

For instance, I am setting as input  a $$ {\tt Reference} (t) = A \sin(\omega t)$$ where $ A = 0.4$ and $\omega = 2\pi\cdot 2KHz$, the outcome of  Discrete-Time
Integrator is my ${\tt Reference} (t)$ integrated a cosine function as expected  $$ve[n] = -B \cos[\Omega n]$$ with $B=2.55*10^4$

However, as shown in the picture  the amplitude $B$ is quite large, which puzzles me. I am aware that the integrator scales the output with $K T_s$, being my $f_s = 200MHz$. So I was wondering whether or not there is an expression that allows me know that scale factor in the amplitude that is happening after the integration, I am aware that is related to my sampling time but I cannot figure out how or what is.

I figured it out coming back to basic maths, I  was making a mistake when integrating my sin function that Matlab block wasn't doing, I was missing behind the $\omega$ variable

$$ve[n] = -B \cdot \cos[\Omega\cdot n] = \int A\sin(\omega\cdot t) = (A/\omega)(-\cos(\omega\cdot t)) + Const$$ $$\mathbf {B = (A/\omega) }$$


Why we take Laplace Transform of functions which converged using Fourier Transform There are several functions for which we know that Fourier Transform will exist but still we calculate its Laplace Transform. Can I know the reason why we need to take Laplace transform for which we know its convergence?

Thanks

There is a large class of functions for which both the Fourier transform and the Laplace transform exist, and for which one can be obtained from the other by setting $s=j\omega$. (Note that even when both exist, the latter need not be the case). So for this class of functions, obtaining the Laplace transform from the Fourier transform (or vice versa) does not require any additional work.

Example: $$\begin{align}x(t)&=e^{-at}u(t),\quad a>0\\
\text{Fourier transform:}\quad X(j\omega)&=\frac{1}{j\omega +a}\\
\text{Laplace transform:}\quad X(s)&=\frac{1}{s +a},\quad |s|> -a\end{align}$$

What the Laplace transform offers is a description in the complex $s$-plane, such as poles and zeros of transfer functions, from which many system properties can be easily deduced. The Fourier transform only shows the behavior on the frequency axis (i.e., the spectrum).


Anisotropic diffusion I am working with speckle noise reduction in ultrasound images. At present I am learning about anisotropic diffusion.I will explain first what I have  understood about anisotropic diffusion.
I learnt the difference between isotropic and anisotropic diffusion using the diagram below:

In case of isotropic diffusion diffusion of the pixel value takes place all across the image.This result in blurring of edges whereas in case of anisotropic diffusion smoothing or diffusion is carried out depending on the image edges and their direction.
It smoothens homogeneous image regions but retains image edges. (I know that the tool that is used to find the strength of the edges and the direction is gradient but I need some more explanation on what is gradient? How we calculate gradient? Please explain me in a lucid manner I have zero knowledge about gradient. Please explain me with easy explanation)
Based on this gradient value we perform smoothing. If the magnitude of the gradient is small then smoothing is done. If the magnitude of the gradient is large then smoothing is stopped.

This graph illustrates the relationship between anisotropic diffusion and gradient magnitude.In the figure the orginal signal is noisy but 2 major peaks are still visible.Applying the diffusison process to the original signal and using the gradient magnitude to attenuate the diffusion process at places where the diffusion process at places
where the edges is strong produce a signal in which the SNR is high.The diffusion at the edges of the signals is attenuated by using a function of the gradient magnitude.(This gradient magnitude equation is what I am unable to understand please see the following equation and explain me in detail. I don't Know why they use divergence and gradient )

( I couldn't able to understand these equations. please explain the equations in relation with the theory above)

Let's proceed step-by-step:

  I need some more explanation on what is gradient? How we calculate gradient? please explain me in a lucid manner I have zero knowledge about gradient.please explain me with easy explanation

a gradient relate to a measure of the derivate, that is, the coefficient of variation. If you think of a measure of position, the derivate is the speed. Now, we do not have a single point, but a field of pixel with values. Now think of it as a landscape, the luminance as the height and the gradient as the slope. It has both a magnitude (how steep it is) and a direction.

In the formulation that you cite, the slope is computed thanks to the divergence operator. In this case it equals to the sum of slope in the directions of the reference frame (that is, in x and y). It is used to allow a more general formulation in non-euclidean spaces.

  This gradient magnitude equation is what I am unable to understand please see the following equation and explain me in detail. I don't Know why they use divergence and gradient

the basic idea behind this equation is to control diffusion by the contours:

diffusion helps to smooth out data by averaging locally. one drawback is that if you do it too long, you lose all information.
contours is the abstract concept of a boundary. Think of it for instance as the contour of a object in an image or the transient in the square wave that you show. To detect the position of these contours, you may use the gradient, and specifically the amplitude of the gradient (how steep the slope is)
now comes the heuristic of anisotropic diffusion: if you diffuse, but that you control the speed of this diffusion by the borders, then, you obtain a better result.

These equations just lay down these concepts quantitatively.


Gabor Filter Parameter I am doing research on images using gabor filter.
I have read the following post:

Apply a Gabor filter to an input image

There must be some formula to determine the value for gamma,sigma,lambda,psi...can i know how to determine it?

The best values for these parameters depend on your application. You should probably try different values on a validation set to find out which ones work best for you.


wavelet analysis/ dimension reduction How wavelet analysis works as a dimension reduction techniques?
The approximations coefficients at higher level of decomposition are the fewer retained coefficients from the originial ?Is it correct?

Well, not really. The wavelet transform (the discrete to be precise) ist just a base transformation of a signal. However, it can be used to compress signals (dimension reduction in some sense), due to the fact that the energy of most natural signals is compactified under the wavelet transform.
The continuous wavelet transform does even the opposite. It blows up your signal and introduces a lot of redundancy.


Texture Classification via DCT How viable would it be to classify the texture of an image using features from a discrete cosine transform? Googling "texture classification dct" only finds a single academic paper on this topic, using a neural network.

For my application, I have a large corpus of labeled images, in which the entire image is a consistent texture (e.g. close-up shots of a blanket, tree bark, a grassy field, etc).

Inspired by a response to a previous question, I was considering the following approach:

split up each image into NxN blocks of pixels
take the DCT of each block
flatten each DCT into a 1xM array and feed it to a K-Means clustering algorithm, and get the cluster label for each DCT
calculate a histogram of clustering labels for each image by counting each label per image from #3
train a SVM classifier by feeding it a set of [(histogram,image label)]

How well would this work? I implemented a similar system, using features extracted via the SIFT/SURF algorithms, but I was only able to get about 60% accuracy.

In what other ways could I use the DCT to classify textures?

So far what you are proposing sounds like a reasonable approach. However, I don't think you will know how well it works until you try it, just like you have tried SIFT.

I have a question though.  Why are you restricting yourself to DCT?  There are lots of representations that have been used for texture classification: co-occurrence matrices, local binary patterns, etc.  The fact that you have only found one paper on using DCT for texture classification would suggest that this is not the most commonly used feature for this problem.  I would recommend that you widen your literature search to see what other features people have used, and how well they have worked.


Does a 20 kHz sine wave sampled at 44.1 kHz need to be reconstituted before playback? I am interested in building a synthesizer but I am confused about the Nyquist-Shannon Theorem and sampling rates.

If I sample a 20 kHz sine wave at 200 kHz and plot the data points, the plot looks similar to that of the true signal. But, many websites I have read reference the Nyquist-Shannon theorem and say audio only needs to be recorded at twice the frequency or higher of the highest frequency one expects to record (or 40+ kHz since humans cannot hear audio frequencies greater than 20 kHz). When I plot the data from sampling the wave at 44.1 kHz, it looks nothing like a sine wave.

Before playing a 20 KHz sine wave sampled at 44.1 kHz on my speakers, must I do something to to reconstitute the original signal from the samples?

What confuses me is that I have seen microcontrollers play audio (.wav files) without reconstructing a signal, simply by taking sample data and playing it at the sample rate -- but I don't think that method will work if the frequency of an audio signal is close to half of the sampling frequency. I have found a lot of literature about sampling rates but nothing about playback rates.

Reminder about Nyquist-Shannon

By the theory of sampling, sampling a signal in the time domain results in periodisation of this signal in the frequency domain, at the period Fs (where Fs is your sampling frequency). Shannon-Nyquist's theorem explains that in order to not loose information when sampling a signal, there shouldn't be overlapping of the periodic repetitions of the signal's spectrum. Considering that the period of repetition is Fs, the total width of your signal's spectrum should exceed Fs.

Moreover, analog signals are real which means their spectrum is symetric relatively to the f=0Hz axis. Thus, if your signal's FT must have a maximum width of Fs, it means that your maximum frequency to avoid overlapping should be Fs/2.

Reconstruction filter

Now that your signal is correctly sampled i.e. there is no overlapping of the spectrum patterns in the frequency domain, you want to reconstruct it in the analog domain for example for audio playback. As explained before, sampling results in a periodisation of the analog signal's spectrum. That means that to go back to a continuous signal, you have to extract only one of the periods of your signal. To do that without distorting your signal, you need an ideal brick-wall filter (infinitely sharp transition band) called a reconstruction filter.
Unfortunately, such a perfect filter is impossible to obtain. Indeed, its impulse response is a cardinal sine (Fourier transform of a gate) which by definition is infinite on both sides of the central lobe.

It can however be approached with very high order FIR filters or convolution with a windowed sinc (this is a pretty deep subject, there's plenty of resources on the internet).

Answer to the original question

Depending on the quality of your reconstruction filter, sampling your signal close to Nyquist will indeed result in a distorted signal. That doesn't mean Nyquist theory is wrong, it is only theory. Depending on your DAC, you often want to oversample your signal if you need very high fidelity. Again, sampling you signal at a higher frequency will result in a periodisation of the signal at a longer period, which basically gives more space for your filter's transition band. The wider the transition band, the more oversampling you will need. As an example, the linear interpolation that your plotting software uses has the frequency shape of a squared sinc. It will thus completely distort your signal.


Autocorrelation of a uniform random process i am currently learning the basics of signal processing.

As you may know the definition of the autocorrelation is different if you look at a random process or for example a deterministic signal

My question is about the autocorrelation of random processes:

Suppose $X$ is a random variable with uniform distribution over $[0,1]$

so: $ f_X(x)=1$ for $(0<x<=1)$

The autocorrelation $r_{XX}(n_1,n_2)$ is defined as: $r_{XX}(n_1,n_2)=E[X(n_1)X(n_2)]$.

if $X$ is stationary up to the second order the autocorrelation is only a function of $\tau$: $r_{XX}(\tau)=E[X(n+\tau)X(n)]$.

If i generate such a random variable in matlab with the "rand" command and compute the autocorrelation (which should be possible because the random process is ergodic [ time and ensemble averages are equal]) i get a strange result which looks more like the convolution of the propability density functions. If i subtract the mean i get the result i would expect, because i assume that X is uncorrelated white noise, so $r_{XX}(\tau)=\sigma_x^2\delta(\tau)$

x=rand(1,100)
Rxx=xcorr(x,x);
subplot(2,1,1)
plot(Rxx);
grid;
title('Autocorrelation function of (X)');

ylabel('Autocorrelation');

y=rand(1,100)-0.5
Ryy=xcorr(y,y);
subplot(2,1,2)
plot(Ryy);
grid;
title('Autocorrelation function of (Y)=(X-0.5) ');

ylabel('Autocorrelation');

So my questions are:

$(1)$ Am i wrong with the assumption that $r_{XX}(\tau)=\sigma_x^2\delta(\tau)$? And if yes, how can we compute $r_{XX}(\tau)$ because we don't know the join-propability densitiy function $f_{X_1X_2}(x_1,x_2)$

$(2)$ Why does the result of the first computed autocorrelation look so strange? (first plot)

As suggested I am adding my comments as an answer.
The sample autocorrelation function (ACF) for $n$ observations is given by

$\hat{p}_x(m) = \frac{\sum_{n=0}^{N-m-1}(x_n - \bar{x})(x_{n+m} - \bar{x})}{\sum_{n=0}^{N-m-1}(x_n - \bar{x})^2} $

To understand why you have a difference between figure 1 and figure 2, lets assume we are looking at observations $x_n$ being zero mean. Let be $y_n = x_n +c$ with c being a constant. Looking now at matlabs help for xcorr for the input $y_n$ we get

$r(m)=\sum_{n=0}^{N-m-1}y_n y_{n+m}^* = \sum_{n=0}^{N-m-1}(x_n + c)(x_{n+m} + c)$

which can be summarized to

$\sum_{n=0}^{N-m-1}(x_nx_{n+m}) + (N-m)c^2$

In the case of xcorr $m$ equals $m=-N,...,N$. Which means you are adding a triangle to your result.

To calculate the ACF matlab has the autocorr function, which takes care of normalization and demeaning of the input.


Applying a 2D band pass filter (Butterworth ?) to an array with NaNs I am trying to filter out the low frequencies of an array looking like this :

For now, I am taking transects of this array and applying a 1D - butterworth band pass filter :

import numpy as np
from scipy.signal import find_peaks, butter, sosfiltfilt, sosfreqz

def Transect_angle(A, angle, Sp, length):
    x0, y0 = Sp[0], Sp[1]
    x1, y1  = x0 + cosd(angle)*length, y0 + sind(angle)*length
    return Array_transect(A, [y0,x0], [y1,x1], 'nearest'), [x0,y0], [x1,y1]

def butter_bandpass(lowcut, highcut, fs, order=5):
        nyq = 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq
        sos = butter(order, [low, high], analog=False, btype='band', output='sos')
        return sos

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        y = sosfiltfilt(sos, data)
        return y

transect, p0, p1 = Transect_angle(Array, angle, start ,length)
tr_withoutnan = transect[~np.isnan(transect)] ### removing nan
filtered = butter_bandpass_filter(tr_withoutnan, low_freq, High_freq, 1)

I would like to apply the same kind of filtering but directly in 2D to my array.
I have two main problems :

For me, filtering is basically a convolution between data and a filter. Based on this, I can directly make the convolution of my array with a butterworth filter having an axial symmetry. However, in my function butter_bandpass_filter, I am using sosfiltfilt, which mention that the filter is appplied forward and backward. I don't know what that means it terms of convolution, but I remember that someaone mentioned it was important in a topic I read a long time ago.
Ho should I deal with the Nans values (in white in my image) ?

Finally, I use a butterworth bandpass filter because I read that it was the filter that was the least likely to produce artefacts. I am open to any suggestion !

filtfilt() is a technique to achieve zero-phase filtering by applying the same filter twice to the data; with the output of the first stage reversed and filtered again in the second stage. Zero phase filtering is a desired property in image processing.
NaN means "not a number" and indicates those indeterminate conditions like $0/0$, $\infty/\infty$, $\infty - \infty$. NaNs cannot be processed arithmetically, hence you should remove all the NaN samples with some suitable values before further processing.


Decomposition of $H(z)$ as maximum-phase, minimum-phase The frequency response is:

$$H(z) = 2-7z^{-1}+7z^{-2}-2z^{-3}$$

I see that it has $3$ zeros: $z_{01} = \frac 12$, $z_{02} = 2$, and $z_{03} = 1$; and $3$ poles in: $$z_x = 0$$

Now, I have to write it like:

$$H(z) = H_{\rm min}(z) H_{\rm max}(z) H_{\rm uc}(z)$$

where $H_{\rm min}(z)$ is the minimum-phase frequency response, $H_{\rm max}(z)$ is the maximum-phase frequency response and $H_{\rm uc}(z)$ only has zeros on $\lvert z\rvert=1$.

For $H_{\rm uc}(z)$, I have one zero on $\lvert z\rvert=1$, so:
$$H_{\rm uc}(z) = z-1$$
Is that OK?
For the minimum-phase frequency, I have all the poles/zeros that are inside the unit circle, so:
$$H_{\rm min}(z) = \frac{z-1/2}{z^3}$$
For the maximum-phase frequency, I have all the poles/zeros that are outside the unit circle, so:
$$H_{\rm max}(z) = z-2$$

But I also know that I can find $H_{\rm max}(z)$ as:
$$H_{\rm max}(z) = H_{\rm min}\left(z^{-1}\right) z^{-M_i}$$
where $M_i$ is the quantity of zeros of $H_{\rm min}(z)$. So:
$$H_{\rm max}(z) = \frac{z^{-1} - 1/2}{z^{-3}} = z^3 \left(z^{-1} - 1/2\right) = z^2 \left(z-1/2\right)$$

So I have $2$ different expressiones for the same $H_{\rm max}(z)$. What am I doing wrong?

So let's check that:
$$H(z) = H_{\rm min}(z) H_{\rm max}(z) H_{\rm uc}(z) \tag{1}$$
where:
$$
H_{\rm min}(z) = \frac{z-1/2}{z^3}\tag{2}
$$
and
$$H_{\rm max}(z) = z-2\tag{3}$$
and
$$H_{\rm uc}(z) = z-1\tag{4}$$

Substituting $(2)$, $(3)$, and $(4)$ into $(1)$:

\begin{align}
H_{\rm min}(z) H_{\rm max}(z) H_{\rm uc}(z) &= \frac{z-1/2}{z^3} (z-2) (z-1)\\
 &= \frac{z-1/2}{z^3} \left(z^2 - 3z + 2\right)\\
&= \frac{1}{z^3} \left(z^3 - 3z^2 + 2z - \frac{1}{2}\left(z^2 - 3z + 2\right)\right) \\
&= \frac{1}{z^3} \left( z^3 - \frac{7}{2} z^2 + 5 z - 1\right)  \\
&\not= H(z)
\end{align}

So something is wrong!

Try:
$$
H_{\rm min}(z) = 1-1/2z^{-1}\tag{2A}
$$
and
$$H_{\rm max}(z) = 1-2z^{-1}\tag{3A}$$
and
$$H_{\rm uc}(z) = 2(1- z^{-1})\tag{4A}$$

Now:
\begin{align}
H_{\rm min}(z) H_{\rm max}(z) H_{\rm uc}(z) &= \left(1-1/2z^{-1}\right) \left(1-2z^{-1}\right) \left(1-z^{-1}\right) 2\\
&= \left(1 - \frac{5}{2} z^{-1} + z^{-2}\right) \left(1-z^{-1}\right) 2\\
&= 2 - 5 z^{-1} + 2z^{-2} - \left(2z^{-1} - 5 z^{-2} + 2 z^{-3}\right)\\
&= 2 - 7 z^{-1} + 7 z^{-2} - 2z^{-3}\\
&= H(z)
\end{align}

I believe your interpretation of Oppenheim and Schafer (OS) is incorrect. And I can find no reference that says the poles have to be outside the unit circle for a maximum phase system. This one and the one I referenced in the comments both only mention the zero locations.

And this is what my copy of OS says:

which also only mentions zero locations.

I believe exercise 5.63 is in error.  The definition used in the body of the book is:

which again does not mention pole locations.

Also, having
$$
H_{\rm max}(z) = H_{\rm min}\left(z^{-1}\right) z^{-M_i}
$$
ensures that $H_{\rm max}$ is causal for FIR $H_{\rm min}$.

So you have two options: either the exercise 5.69 definition is wrong, or the relationship between $H_{\rm max}$ and $H_{\rm min}$ is wrong.


What does convolution has the meaning? As I know, if we want to know the LTI system output, then we do convolution between input x[n] and impulse response h[n].
but actually,in this question, I want to know what does convolution has the behind meaning.

Why do we do convolution(sum of products) not just using adder or multiplier for calculation between input signal and impulse response?

I explained in "Why is time inversion?"

Red is a reference delta impulse (of height 1), whereas green is a typical response to that impulse, denoted by function h(t). In LTI, output is proportional to the input, that is if we have a delta impulses at the input at time 0, output at time t will be a*h(t). Now, instead of single impulse at the origin, you apply a series of input impulses at various times. What will be the output? Say, there was input impulse of height $a_1$ at $T_1$. Since current time is $t$, impulse occurred $t-T_1$ seconds ago and its contribution to the current output y(t) is $a_1 h(t - T_1)$. There is another contribution from another impulse, occurred at $T_2$. Its contribution is $a_2 h(t-T_2)$. So, $y(t) = a_1 h(t-T_1) + a_2 h(t-T_2)$. You simply add up the contributions beacuse of LTI linearity.

In general, you have $y_t = \sum_0^t {a_i h(t-i)}$. That is a convolution formula.

It also appears when you multiply to polynomials $(a_0 + a_1 z + a_2 z^2 + ...)(b_0 + b_1 z + b_2 z^2 + ...) =\sum_0^\infty {c_n z^n} $ where $c_n = a_i b_{n-i}$. That is why you tend to represent series as z-transforms. In this case you can simply multiply them and have convolution on the background.


How do I process IQ data received from a signal analyzer I used a signal analyzer to get IQ sample data from a blue tooth device. The IQ data received is in the form of a .MAT file and I examined the file contents and see an + ibn samples where a, b are constants, n is the sample number and the format is a complex number. I was told the real part is the I and the imaginary part is the Q. Are these samples received baseband I Q data? I am trying to determine the carrier frequency offset. I was trying to determine the frequency of the signal. I felt I could measure the phase difference between samples per time. However, I have the angles being only between pi/2 and -pi/2 when I do arctan(b/a). I assume I should unwrap the signal phase. In MATLAB there is the unwrap command (var = unwrap(phase)). But that does nothing to the values I received. I am unsure of how to proceed with my search for the carrier frequency offset. Any help?

If you're looking for the phase progression I would recommend not using atan() in MATLAB to extract the phase angle.  Use the built-in function angle() to extract the phase from the complex samples properly accounting for the quadrant of the complex plane in which each sample appears.  This resulting phase can then be unwrapped.

unwrappedPhaseRad = unwrap(angle(inputComplexSamples));

where unwrappedPhaseRad is the unwrapped phase in radians and inputComplexSamples is a vector (1 dimensional array) of your complex input data samples.


Channel Shortening for underwater acoustic channel communication I'm new in this forum, hope that we can learn from each other well.  I need to ask about the channel shortening for underwater acoustic channel.

what are the common algorithms in that field? is there a common source code in MATLAB which is used to be compared with such new algorithms?

I don't see many papers in that field, why?

please help me,

thank you so much.

When I see "channel shortening," I think of DFE coefficient calculations, regardless of the application. This is used in cases from simple DFE equalization to channel shortening for xDSL with DMT to lessen CP requirements or for RSSE in single-carrier systems. The point is not to actually do a DFE; rather, the point is to use the coefficient calculation and paradigm.

The classical DFE architecture is to have two filters: a feedforward (FF) filter and feedback (FB) filter with a slicer after the FB subtraction and decision-directed feedback passed through the FB filter. The intent is to use a FF filter to compress the residual channel response into a smaller number of taps than you started with. For DFE applications, the intent is also to transform a possibly non-mininum phase response to minimum phase to reduce error propagation. The desired structure is to have a FF filter that is (close to) all-pass and simply transforming phase, while the FB taps of the DFE represent the residual channel. That said, the concept of compressing the channel into a shortened response can generalize to the case where the FF filter is not all-pass. In this case, it looks more like a FF equalizer for a partial-response system.

I suggest investigating DFE coefficient calculation and motivations to see approaches for channel shortening. It may also be instructive to look at partial-response literature. In the case of PR, you may see descriptions of equalization to a predetermined response, but it can also be used more generally for compression of a response into a reduced number of taps.


Matlab: if given nothing but an audio file, how do I modify it to produce incrementally lower signal to noise ratios? I have an audio file with a high sampling rate of 250khz. It contains animal noises. I've built a program which counts all the animal noises.

However, I want to find out how robust the program is and find what the minimum signal to noise is for the program to still work.

Hence the plan was to modify my audio file to produce different (incrementally lower) signal to noise ratios.

I would run the program for all of them and since I know how many the program should count, when I start getting lower counts, I will what the minimum signal to noise ratio is.

however I have no idea how to go about doing this. I still don't even know how to find the signal to noise of my audio as it is. my guess was to assume the background noise is similar throughout and hence, copy a section of the audio which has no animal calls, just noise, then use matlabs snr command:

r=snr(audio,backgroundnoise);

can anyone help with how to program this?

It is not enough to just randomly add noise, I need to know the actual signal to noise values in order to be able to quantify the minimum SNR needed for the program to still work

I like your noise model, I think that is really the best way to find SNR if you only have access to recordings. Using this same model couldn't you take your backgroundnoise scale it (duplicate it or trim it as necessary) and then add it to the signal? In this sense the noise wouldn't be entirely random but somewhat system specific

edit

in order to calculate your new noise value you would simply rearrange the SNR calculation solving for noise. Then scale the original noise accordingly. you can use some code like this

desired_snr_db = 10;    %change this value

est_signal = audio - backgroundnoise;

%the original matlab equation for SNR is:
%signalPow = rssq(x(:)).^2;
%noisePow = rssq(y(:)).^2;
%r = 10 * log10(signalPow / noisePow);
%solving for noisePow we get
noisePow = rssq(est_signal(:)).^2 / (10^(desired_snr_db / 10));

%assuming rssq(new_noise).^2 = (scale_factor * rssq(old)).^2
scale_factor = sqrt(noisePow / rssq(backgroundnoise).^2);

new_noise = scale_factor * backgroundnoise;

%add noise to our estimated signal
new_signal = est_signal + new_noise; %you may need to do stretching
                                     %trimming to get the right sizes

new_snr = snr(new_signal - new_noise, new_noise)

if you want you can double check the snr by finding a lull in the new signal and comparing it to the known sounds just as before


Getting frequencies out of a FFT I use this FFT and wrote a short program to test it.

package test_FFT;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.TransformType;

public class FFT2 {

public static void main(String[] args) {
    double[][] array=new double[2][8];
    //An array of complex numbers with real part array[0][*]
    //and complex part array[1][*]
    //The output of the transformation will be saved in the input array
    array[0][0]=5.0;
    array[0][1]=2.0;
    array[0][2]=3.0;
    array[0][3]=4.0;
    array[0][4]=5.0;
    array[0][5]=6.0;
    array[0][6]=7.0;
    array[0][7]=9.0;
    FastFourierTransformer.transformInPlace(array,DftNormalization.STANDARD,TransformType.FORWARD);

    for(int i=0;i<8;i++){
        System.out.println("real "+array[0][i]);
        System.out.println("complex "+array[1][i]);
    }

}

}

In the output array I get amplitudes of sin and cos functions. The information about the frequencies should depend on the position within the output array.
After some research on this page, I still don't understand how to calculate frequencies out of array positions.

I learned that there are many flavours of how to perform a FFT. Has anyone of you detailed knowledge on how calculate frequencies out of output of the FFT I use? A code sample computing frequencies for the output of the example below would be greatly appreciated. Thank you

At this moment i'm not hos i can't check my answer...
You say

  In the output array I get amplitudes of sin and cos functions.

Well that is you frequency array. Probably not the way want it. I suspect youre looking for the magnitude and maybe the phase.

the magnitude is quite easy.. it is the $abs(array (real + imag)$.
the phase $ \phi = \arctan \frac{imag}{real} $

in this way it is common to show the half the frequency array because of it's symmetry.
if you plot $array[0]$ to $array[N/2]$ you shoult get your frequecy array from $ \frac{1}{N}$ to $\frac{F_s}{2}$


Is there a version of Welch's method that doesn't look for power? Welch's method splits a time signal, $x(n)$ into $M$ periodograms $P_m$,

$P_{x_m,M }(k) = \frac{1}{M}|F_k(x_m)|^2$

and averages them to give the Power Spectral Density (PSD),

$S_{x}(k) = \frac{1}{K}\sum_{m=0}^{K-1} P_{x_m,M }(k)$

But if I am not interested in the power can I use something like an averaged Fourier transform to represent the Fourier transform for the whole signal $x(n)$?

I.e. The exact same as above but instead of averaging the periodograms, $|F_k(x_m)|^2$, I just average the fourier transforms $F_k(x_m)$ within each window?

Kind of like a STFT except averaging the windows, and is there any merit to using something like this?

Since an FFT is a linear operator, adding up the complex results of a sequence of FFTs of short windows is the same as doing a single short FFT on the the vector addition of all those short windows.

Note that for signals that are exactly integer periodic in the FFT width (sequential 0% overlapped windows), the vector addition will constructively interfere.  But for signals that are not exact integer periodic in the FFT width, there can be significant destructive interference.  So the result will be highly scalloped, with a narrow high gain bandpass around each FFT bin center, and potentially a lot of signal loss between bin centers.  Sometimes this is useful, sometimes not, but a very different result from a single longer FFT over the entire window set.


Problems designing high pass FIR filter with $f_\textrm{cutoff} \ll f_\textrm{s}$ Is there any inherant problems with high-passing a signal with a very low cutoff frequency compared to the sampling rate?

I'm not sure what I am missing as I cant find any info on it.

Examples with gives bad results:

sample_rate = 250
nyq_rate = sample_rate / 2.0
f_hp_freq = 0.1
f_hp_transition = f_hp_freq + 0.1
numtaps = 125
b = signal.remez(numtaps, [0, f_hp_freq, f_hp_transition, nyq_rate], [0, 1], Hz=sample_rate, type='bandpass')

or

b = signal.firwin2(125, [0, f_hp_freq, f_hp_transition, nyq_rate], [0, 0, 1, 1], nyq=nyq_rate)

by reducing the sampling rate to less than $10\textrm{ Hz}$ or by significantly increasing the filter order < 1000 I can get acceptable results, but no other way.
What am I overlooking?

@poipoikuroi: I don't think you are "overlooking" anything. If IIR filters are allowed you can build a simple "DC Blocking" filter with a narrow stopband (using floating-point arithmetic) centered at zero Hz.  (Search the web for DC Blocker.) However, those filters have wildly nonlinear phase near zero Hz.

If linear phase is mandatory then your super-narrow stopband width requirement makes your FIR filter impractical due to needing thousands of taps.  So as you've noticed, you'll have to either widen your stopband width or lower your data sample rate. I developed a computationally-efficient linear-phase DC blocking filter that uses cascaded moving average filters. Perhaps it might be of interest to you. See:
https://www.dsprelated.com/showarticle/58.php

Does your eeg software package use one or two FIR filters to remove both DC and line noise? If your eeg software uses a single high-pass FIR to remove both DC and line noise, then I'll bet that filter's stopband width is, say, 90 or 100 Hz. (Which is much more "practical" than the FIR filter you're trying to design.)


STFT Window Size is 2048, Why is the output 1025? I feel like I am having a brainfart over here and can't seem to remember what's going on with STFT outputs.

Consider these two lines of code from the Python library Librosa:

# Window the time series.
y_frames = util.frame(y, frame_length=n_fft, hop_length=hop_length)

# Pre-allocate the STFT matrix
stft_matrix = np.empty((int(1 + n_fft // 2), y_frames.shape[1]),
                       dtype=dtype,
                       order='F')

In the first line, we have a function that creates a matrix with a window length of n_fft (2048).

Then in the next line, we pre-allocate our STFT, but our window length is now 1025 instead of 1024 as dictated by the 1+n_fft // 2?  Where does this extra frequency bin come from?  Why is not just 1024?

The continuous Fourier transform possesses symmetries  when computed on real signals (Hermitian symmetry). The discrete version, an FFT (of even length) possesses a slighty twisted symmetry.

The DC coefficient ($F(0)$) is real, as well as the Nyquist one ($F(N/2)$). In between, you get $\frac{2048-2}{2}=1023$ "complex" coefficients, "duplicated" in positive and negative frequencies.

So for real signal, each STFT frame can be represented by $1023+2$ frequency bins, the remaining   1023 being recovered by Hermitian symmetry.

You can get complementary information at FFT of Pure Real Sequences:

  As a result, we can see that if N is even, both $F(0)$ and $F(N/2)$ must
  be real. Given these two values and the complex values $F(1)...F(N/2-1)$,
  (I.E. N numbers in total) the sequence is completely characterised.


Is ICA appropriate for separating mixed signals when all source signals are NOT detectable at all sensors? A generic implementation of ICA for the separation of a mixture of $N$ signals into their $M$ constituent components requires that the signals be assumed to be a linear instantaneous mixture of the sources.  Every description of ICA that I have come across seems to take for granted the fact that all $M$ sources are present to some extent in all $N$ signal mixtures.

My question is, what if the $M$ sources are only present in some but not all of the signal mixtures?

Does this scenario violate the fundamental assumptions necessary for ICA to be able to separate these signals?  (Assume, for the sake of argument, that we are dealing with an overcomplete or complete system ($N>M$ or $N=M$), and that each of the $M$ source signals are in fact statistically independent from each other).

The implementation I am considering using ICA for, in which this situation arises, is the following: I have data from 4 different types of sensors, each with a different number of channels.  Specifically, I have 24 channels of EEG data, 3 channels of electrooculographic (EOG) data, 4 channels of EMG data, and 1 channel of ECG data.  All data is recorded simultaneously.

I would like to identify the contributions of the ECG, EMG, and EOG signals within the EEG data so that I can remove them.  The expectation is that EMG+ECG+EOG signals will be picked up by the EEG sensors, but not vice versa.  Also, EOG and EMG will likely contaminate each other and be contaminated by ECG, but ECG will probably be pretty isolated from all the other signals.  Also, I'm assuming that where mixing does occur, it is linear and instantaneous.

My intuition tells me that, hypothetically, ICA should be smart enough to return mixing filters with very small (close to 0) coefficients to account for the lack of contribution of sources to a mixed signal.  But I am worried that something about the way that ICA demixes the signals inherently enforces the expectation that all sources will be present in all mixtures.  The implementation I am using is FastICA, which is a projection pursuit-based approach.

You should be fine, zeros in the mixing matrix are not a problem....and theoretically it should converge even quicker than if all sources existed in all sensors.


How can I increase the zero-padding on the DFT? Now I'm trying to understand and implement the zero-padding work in my example.
As I know, meaningless zero values are affect to better DFT result. So I want to see how zero-padding works.

clear all
clc

fs = 50; %50Hz
T = 1/fs;
N = 8;
t = (0:N-1)*T;

y = 5+cos(2*pi*12.5*t)+sin(2*pi*18.75*t);

X=zeros(1,N);

for k = 0:N-1
    for n = 0:N-1
        X(k+1) = X(k+1) + y(n+1)*exp(-j*(2*pi/N)*k*n);
    end
end

x_mag = abs(X);

f_2 = fs*(0:(N-1))/N;
stem(f_2,x_mag,'-ok','linewidth',2,'MarkerFaceColor','black');

For N=8:

For N=80:

But I feel hard to implement the increasing zero padding in my example code.
Would you please show me how I can increase the zero-pad?

I found the better result when I increased the $N$ number to 80 from 8.
But I'm not sure if increasing $N$ means increasing the zero-pad.

Your code does not implement the zero padding method. Below is the corrected version that does. Note that zero padding does not improve frequency resolution. It mereley provides interpolated samples of the original spectrum.

   clear all
clc

fs = 50; %50Hz
T = 1/fs;
N = 8;
t = (0:N-1)*T;

y = 5+cos(2*pi*12.5*t)+sin(2*pi*18.75*t);

M = 8;
K = M*N;       % DFT increased size
yk = [y zeros(1,(M-1)*N) ];  %zero padded input

X=zeros(1,K);   % high density

for k = 0:K-1
    for n = 0:K-1
        X(k+1) = X(k+1) + yk(n+1)*exp(-j*(2*pi/K)*k*n);
    end
end

x_mag = abs(X);

f_2 = fs*(0:(K-1))/K;
stem(f_2,x_mag,'-ok','linewidth',2,'MarkerFaceColor','black');


Create distortion from basic linear (and non-linear if neccessary) DSP elements I'm studying mechatronics and I'm intrested in DSP basics. My lecturer said that there are four basic linear DSP elements:

Adder (and other mathematical operations)
Amplification (shown on diagrams as triangle)
Delay ($z^{-1}$)
Junction

I know from electronics that I can create distortion using for example op-amps and clipping circuits - and I can do it practically. But is there a way to create distortion using four DSP elements theoretically?
I know distortion is non-linear - so if not this four elements, what kind of elements existing in theoretically DSP should i use?

Using those four basic elements will allow you to implement linear systems, which can change the magnitude and phase of the input signal, but which will not add the harmonics that are expected from a distortion effect. In order to create distortion in that sense (i.e., non-linear distortion) you will need some non-linear element. The most basic implementation would use a (soft) clipping function, possibly followed by some linear filter to shape the spectrum.

A very accessible reference is the DAFX (Digital Audio Effects) book edited by Udo Zölzer. It has a chapter on non-linear processing, including distortion and overdrive. Note that I found an error in one of the equations for producing distortion (Eq. (5.9)). You can find the corrected formula in this answer to a related question.


Using the system function $H(z)$ to the solve the output of a DT LTI system I have a question which asks me to find the magnitude and phase response of

$$
H(z) = \frac{1}{1-0.5z^{-1}}
$$
for which I solved by doing
$$
H(e^{j\Omega}) = \frac{1}{1-0.5e^{-j\Omega}} = \frac{1}{1-0.5\cos(\Omega)+j0.5\sin(\Omega)}
$$
$$
|{H(e^{j\Omega})}| = \frac{1}{\sqrt{(1-0.5\cos(\Omega))^{2}+(0.5\sin(\Omega))^{2}}} = \frac{1}{\sqrt{\frac{5}{4} - \cos(\Omega)}}
$$
$$
\angle H(e^{j\Omega}) = -\tan^{-1} \left(\frac{\frac{1}{2}\sin(\Omega)}{1-\frac{1}{2}\cos(\Omega)}\right)
$$
But, I am unsure how to solve the second part of the problem, which asks me to compute the system response to
$$
f[k] = A\cos\left(\frac{\pi}{2}k +\frac{\pi}{10}\right)
$$

Let's first briefly recall the meaning of $H(z)$.
Let $y[n]$ be the output of the LTI system. Given an input $x[n]$, the output can be written
$$y[n] = h[n] *x[n] = \sum_{k=-\infty}^\infty h[k]x[n-k]$$
where $h[n]$ is the impulse response of the LTI system en $*$ represents the convolution. Let's now assume that $x[n]$, the input of the LTI system, is of the form $z^n$. In this case,
$$y[n] = \sum_{k=-\infty}^\infty h[k]z^{n-k} = z^n\left(\sum_{k=-\infty}^\infty h[k]z^{-k}\right) = z^nH(z)$$
where $H(z)$ denotes the transfer function of the LTI system (the z-transform of the impulse response). It can be seen that that $z^n = e^{j\Omega}$ is an eigenfunction of the system, that is, it is only affected by a complex factor (i.e. $H(z)$, the so-called eigenvalue) at the output of your system.

Now, your input $f[n]$ can be written as
$$f[n] = \mathcal{R}\{Ae^{j(\frac{\pi}{2}n + \frac{2\pi}{10})}\} = \mathcal{R}\{Ae^{j\frac{2\pi}{10}}e^{j\frac{\pi}{2}n}\}$$
and $Ae^{j\frac{2\pi}{10}}e^{j\frac{\pi}{2}n}$ matches the form of an eigenfunction (with a complex factor $Ae^{j\frac{2\pi}{10}}$). The output of your system to $Ae^{j\frac{2\pi}{10}}e^{j\frac{\pi}{2}n}$ is then given, using linearity of your system and the polar notation of $H$, by
$$y[n] = Ae^{j\frac{2\pi}{10}}H(e^{j\frac{\pi}{2}})e^{j\frac{\pi}{2}n} = Ae^{j\frac{2\pi}{10}}|H(e^{j\frac{\pi}{2}})|e^{\arg\{H(e^{j\frac{\pi}{2}})\}}e^{j\frac{\pi}{2}n}.$$
Finally, again by linearity, the output of your system to $f[n]$ is given by the real part of that.


OpenCV/C++ connect nearby contours based on distance between them I have to connect nearby contours in an image based on the distance between them which specifies if the contours are to be connected.

Now there's a question already on the same problem here https://stackoverflow.com/questions/8973017/opencv-c-obj-c-connect-nearby-contours but here he merges all contours into a single one. This I don't want. I don't think that there is some function in opencv for this but you can suggest an algorithm for that. My application goes like this:

I'm detecting hands, so I used a skin detection algorithm to determine them but since my skin ain't white and maybe because of lightening conditions sometime the contour breaks of at elbow. So i want the nearby contours to be connected but not all of them (because both both my hands will be there in contours.) (By hands I mean from shoulder to palm.)

Moreover I think that by using some edge detection I will get my hands boundaries and detect if some of this patch inside this boundary is detected as skin then whole region within this boundary will be detected as skin but i'm not sure how to do this part.

Any Help will be appreciated. Thanks in advance

Sample image:

In this Image I wanna connect points (8 connectivity) which are less than say 40 pixels of distance so that I'll get my left hand as a single contour

My aim is to only get hand's contour (I don't care about any other region)

If you are not worried about the speed or exact contour of hand, below is a simple solution.

The method is like this : You take each contour and find distance to other contours. If distance is less than 50, they are nearby and you put them together. If not, they are put as different.

So checking distance to each contour is a time consuming process. Takes a few seconds. So no way you can do it real time.

Also, to join contours, I put them in a single set and drew a convex hull for that set. So the result you are getting is actually a convex hull of hand, not real hand.

Below is my piece of code in OpenCV-Python. I haven't gone for any optimization, just wanted it to work, that's all. If it solves your problem, go for optimization.

import cv2
import numpy as np

def find_if_close(cnt1,cnt2):
    row1,row2 = cnt1.shape[0],cnt2.shape[0]
    for i in xrange(row1):
        for j in xrange(row2):
            dist = np.linalg.norm(cnt1[i]-cnt2[j])
            if abs(dist) < 50 :
                return True
            elif i==row1-1 and j==row2-1:
                return False

img = cv2.imread('dspcnt.jpg')
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
ret,thresh = cv2.threshold(gray,127,255,0)
contours,hier = cv2.findContours(thresh,cv2.RETR_EXTERNAL,2)

LENGTH = len(contours)
status = np.zeros((LENGTH,1))

for i,cnt1 in enumerate(contours):
    x = i
    if i != LENGTH-1:
        for j,cnt2 in enumerate(contours[i+1:]):
            x = x+1
            dist = find_if_close(cnt1,cnt2)
            if dist == True:
                val = min(status[i],status[x])
                status[x] = status[i] = val
            else:
                if status[x]==status[i]:
                    status[x] = i+1

unified = []
maximum = int(status.max())+1
for i in xrange(maximum):
    pos = np.where(status==i)[0]
    if pos.size != 0:
        cont = np.vstack(contours[i] for i in pos)
        hull = cv2.convexHull(cont)
        unified.append(hull)

cv2.drawContours(img,unified,-1,(0,255,0),2)
cv2.drawContours(thresh,unified,-1,255,-1)

Below are the results i got :


Plotting wavelet co-efficients I need to take the wavelet transform of an image (for reference, I will be using Debauchies wavelets, and the test image will either be Lenna or matlabs default woman image),and plot the coefficients in order of their magnitude. This is as an example for a document on compressive sensing I'm writing.

I'm a little unsure of of which wavelet function to use, and then what to plot to best show that many of the coefficients are zero. So far I've got:

load woman
[c,s] = wavedec2(X,2,'db1');
a = sort(c)
plot(a)

Which does have a nice power law distribution:

However, I'm still a little confused and unsure if this image shows what I think it shows.

Could someone suggest an alternative?

You are on the right track but need perhaps to quantify your goal in mathematical terms.

The sorted coefficients is incidently also the inverse of the cumulative function (see for instance the eq 6 of this paper).

Then to quantify the sparsity of your signal you could

compute some measure using the cumulative distribution function, such as the ratio of coefficients above a threshold (the 5% max quantile),
and compare with a control condition, like doing the same process
with a random signal


Impulse response time period in circular convolution While considering an input to be periodic of Period N, can the impulse response not be periodic of period greater than N ? If it can be, how can one compute it’s convolution?

Circular convolution assumes that all signals ($x[n]$, $h[n]$ and $y[n]$) are periodic in the same integer $L$. When any of the signals are shorther than $L$ then they are padded with enough zeros to make them periodic with $L$. When $x$ or $h$ has a period larger than $L$ then there will be aliasing in the computed output $y[n]$ which is still periodic with $L$.


Applying Kalman filter to a data set I went through the answer Kalman filter in practice and it seems we must know all the first and second order properties of random variables to apply the Kalman filter.

But when I only have a set of data $y_1,\dots,y_n$ and a given model $G$ and $F$ with assumed noise processes, could someone please explain to me how to apply Kalman filter?

Thanks.

Signal Model

The first thing you really need when you start to use a Kalman filter is what I would call the "signal model".  That means, you need to take a guess (and sometimes that's all it is, though usually you can have an educated guess) at how the signal you're measuring, $y$, was generated.

The usual signal model comprises $F$, $G$, $H$, $Q$ and $R$:

$$
x_{k+1} = F x_k + G u_k + w_k\\
y_k = H x_k + v_k
$$

where $Q$ is the covariance of $w_k$ and $R$ is the covariance of $v_k$.

What are you trying to achieve?

Implicit in the signal model is that you've chose the "state variables" (the components of $x_k$) that you are working with.  Usually, what you are trying to achieve is to measure something and, from that, infer the value of one or more of the state variables.

What's missing?

Based on what you've said, the only thing that seems to be missing is how you get from the states $x_k$ to your measurements $y$ (i.e. $H$).


Analyzing a particular discrete-time LTI system for input signal $x[n]=(1/3)^n$ for *all* $n$ I'm considering the following problem from some course notes.

  Suppose the following is known about a discrete-time LTI system:

  Given the input $x[n]=(1/3)^n$ for all $n$, the system produces output $y[n]=(2/3)(1/3)^n$ for all $n$.
  Given the input $x[n]=(1/2)^nu[n]$, where $u[n]$ is the unit step, the output is $y[n]=\delta[n]+a(1/4)^nu[n]$.

  The problem asks for the value of $a$ and for the response $y[n]$
  given the input $x[n]=(1/6)^nu[n]$.

I thought from (1) I could determine the system function by taking a Fourier or Z transform, and then use that with (2) to deduce the value of $a$. For example, the Z transform of (2) says $Y(z)=H(z)X(z)$, where $Y$ and $X$ are easily computed (or obtainable from a table). If I could deduce $H$ from (1), I could get $a$.

But, as discussed here, the signals $(1/3)^n$ and $(2/3)(1/3)^n$ for all $n$ do not have Fourier or Z transforms, so I'm stuck. I could imagine there being a typo: perhaps $(1/3)^n$ should be $(1/3)^{|n|}$. But if that's not the case, I'm not sure how to proceed.

This problem tests the understanding of two basic properties of (discrete-time) LTI systems. First, for an input signal of the form $x[n]=z_0^n$, $-\infty<n<\infty$, the output is given by

$$y[n]=z_0^nH(z_0)\tag{1}$$

where $H(z)$ is the system's transfer function, and where it is assumed that $z_0$ is inside the region of convergence (ROC) of $H(z)$. [Note that $(1)$ is not true for $z_0$ outside the ROC.]

The second property is simply the $\mathcal{Z}$-domain input-output relationship

$$Y(z)=X(z)H(z)\tag{2}$$

Applying $(1)$ to your example you can conclude that

$$H\left(\frac{1}{3}\right)=\frac23\tag{3}$$

and from $(2)$ you can get an expression for $H(z)$ in terms of the constant $a$. This constant can be determined from Eq. $(3)$.


Convolution in the frequency domain We often hear that "convolution in time is the same as multiplication in frequency", and vice versa, that "convolution in frequency is the same as multiplication in time".

So in a typical windowing operation, we do a point-wise multiplication of a signal $x[n]$, with a window $w[n]$. This means that in the frequency domain, we are performing $X(f) *
W(f)$.

My question is the following: I wish to illustrate this for myself, but I am not sure what lengths FFTs to take, and what type of convolution to do. (Circular? Linear? etc). Those details seem to be missing from all sources I look at.

So, if you want to do convolution in the frequency domain, what lengths FFT do I pick to compute $X(f)$ and $W(f)$, and what type of convolution am I doing? Thank you.

the FFT is a fast method of doing the DFT.    yes, multiplication in one domain corresponds to convolution in the reciprocal domain.  and since the DFT is a totally circular operation, any consequential convolution done in one domain by the DFT by multiplication in the other domain is circular convolution.  you don't have a choice about it.  to make a circular convolution look like linear convolution, you must zero-pad sufficiently.

how many zeros to pad?  need to consider how fast convolution (overlap-add or overlap-save) works.  if $N$ is the DFT length, the non-zero lengths of the two sequences being convolved, must not add to more than $N+1$ otherwise there will be overlapping and aliasing of non-zero components.


Convolution of two large signals in MATLAB I need help regarding convolution. I am trying to convolve a signal which is $15\textrm{ ms}$ long, and the sampling frequency is $10\textrm{ GHz}$. I need to convolve it with a rectangular pulse, which has $15\textrm{$\ \mu$s}$ pulse width. When I am trying to use conv function in MATLAB, my PC goes busy state, and it does not give me any answer. Can anybody help me with this problem, please?

If your shorther pulse signal is really a rectengular waveform then please look for the other answer, but otherwise for a general waveform pulse the following code snipped (excerpted from Maximillian's previous post) shows the actual results (on my laptop with MATLAB R2015) of timings for a frequency domain vs time domain implementation of the convolution operation in simplest terms.

Note that I have reduced the sampling rate Fs from 10 Ghz to 1Ghz due to memory reasons and also I have slightly adjusted the sequence lengths to result in a power of $2$ FFT length so that FFT could run at full speed. (this is slightly unfair for a practical course but here I would like to see the maximum attainable speed gain, hence it would be ill posed to use inefficient FFT lengths for such a purpose)

Fs = 10E8; % reduced for MEMORY reasons
T = 17e-3; % slightly adjusted
t = 0:1/Fs:T;
t = t(1:(2^24-2^14+1)); % its length is adjusted to create a 2^N factoring later.
y = sin(2*pi*1000*t); % The input signal
L1 = length(y);

w_rect = 17e-6;  % width of rect
R = ones(1,Fs*w_rect);
R = pulse(1:2^14);
L2 = length(R);

tic
y2 = real(ifft(fft(y,L1+L2-1).*fft(R,L1+L2-1)));
tdft = toc

tic
y3 = conv(y,pulse);
tconv = toc

eff = tconv/tdft

The result is the following: tdft = 1.4769 s, tconv = 42.14 s and hence speed gain = eff = 28.5x.

On the other hand, for slightly different but otherwise unfactorable FFT lengths, the speed gains decreased down to 10x as understandable.

This is exactly the situation where you should seriously consider implementing your convolution in the frequency domain (using efficient FFTs), rather than in the time domain!

The condition is triggered by the fact that the shorter pulse signal (with about 150k samples) is already too long for a direct time domain convolution.

Lets make a brief comparison between the two approaches for the signals of lengths about $L_1 = 150 \times 10^6$ samples and $L_2 = 150 \times 10^3$ samples:

In the first case a direct time domain implementation would approximately require $L_1 \times L_2$ many real MACS which is about $(150 \times 10^6) \times (150 \times 10^3) = 2.25 \times 10^{13}$ real MACS.

Now compare this to an FFT based frequency domain implementation (according to the convolution in time domain equals multiplication in frequency domain theorem of DFT). Again using rough numbers there will be three steps: 1- Convert both signals to frequency domain of the same (and sufficiently long) lengths (which is actually $L_1+L_2-1$ )using efficient FFT, 2-Multpily the obtained DFTs and get the output DFT, 3-Convert the result back to the time domain using the inverse FFT.

The first step requires about $2 \times (150 \times 10^6) \times \log_2(150*10^6) \times 4$ real MACS which is about $3.25 \times 10^{10}$ real MACs

The second step would require roughly $4 \times (150 \times 10^6) = 6 \times 10^8$ real MACs

And the final, third, step would require about $1.6 \times 10^{10}$ real MACS.

Hence the totality of the FFT based frequency domain convolution implementation would require roughly: $(3.2 \times 10^{10}) + (6 \times 10^8) + (1.6 \times 10^{10}) = 5 \times 10^{10}$ real MACS.

Now compare this to the time domain implementation which required $2.25 \times 10^{13}$ real MACS. This could potentially result in about 500x speed improvement assuming a basic implementation. So for example if the convolution would take 5 hours in direct time domain implementation then it should take about 36 seconds in DFT based case!

So you better follow the second approach. That being said, however, now current available multicore CPU + GPU based extensively vectorized (parallel) implementations can dramatically reduce the direct time domain computation time of convolutions, however again then the efficiency of FFT will also benefit from this parallelization and therefore you should, nevertheless, see an order of magnitude improvement in speed when using a DFT approach.

PLease also consider that such long sequences can be processed in blocks and that also would account for certain imprevements in memory and speed performances.

EDIT: I've tried to implement the DFT processing in an non-blocked based mode, where the whole signal processed at once in memory but that required too much RAM (about at least 10 GBs for MATLAB alone in my PC which has only 8GBs of RAM for all programs currently so it started paging and went into freezing, (the computer simply stcuk without any response even to ctrl+alt+del, so I had to hard reset it) If you have more memory such as 16 GB or above you may still try and it may probably do ok, as I assumed it would at most require about 10 GBs of RAM (plus any other programs using your remainig RAM). But otherwise you should consider dividing your larger signal (which 150 M samples) into smaller pieces or use some other technique.


Convolving time reversed discrete input with time reversed discrete impulse response Question 3.1 d from Chapter 3 of Oppenheim's Signals & Systems
I have

$$x[n]=(-1)^n{u[-n]-u[-n-8]}$$
$$h[n]=u[n]-u[n-8]$$

and the question asks for $y[n]$ which is the convolution of $x[n]$ and $h[n]$ namely

$$y[n]=x[n]\star h[n]$$

I'm stuck with the limits on the sum

$$\sum_{k=-\infty}^{\infty}{\left((-1)^{k}(u[-k]-u[-k-8])\right)\left(u[n-k]-u[(n-8)-k]\right)}$$

Should the limit be '$0$' or '$n-k$' or '$n-8$'?

Here is another example of how incredibly poor choice of notation misleads students everywhere. I do not possess a copy of any of the various tomes on signal processing that seem to be revered as the fifth Gospel on this site, and so
I do not know if the book cited by the OP actually has written $(1)$,
but I assert that it is against mathematical common sense to write things like
$$y[n] = x[n] \star h[n].\tag{1}$$
It is generally accepted mathematical convention that is a symbol
has the same meaning everywhere it appears in an equation or expression.
For example, when we write
$$y(t) = x(t)h(t),~ -\infty < t < \infty, \tag{2}$$
we mean by this that for every choice of real number $t$, the value
of $y(t)$ is the same as the product of the values of $x(t)$ and $h(t)$;
for example, $y(3) = x(3)h(3)$, and $y(313.012) = x(313.012)h(313.012)$,
ans so on and so forth.
So how are we supposed to make sense of this monstrosity $(1)$?
Clearly $y[3]$ is not $x[3]\star h[3]$ whatever meaning we might
ascribe to the latter quantity; the value of $y[3]$ depends on
a lot of other values of $x$ and $h$.  The poor notation misleads
students into writing things like

$$\sum_{k=-\infty}^{\infty}{\left((-1)^{k}(u[-k]-u[-k-8])\right)\left(u[-n-k]-u[(-n-8)-k]\right)}$$

which is not the correct expression for the convolution (or the correlation,
for that matter) of
$$x[n]=(-1)^n{u[-n]-u[-n-8]}\quad \text{and} \quad
h[n]=u[n]-u[n-8]\tag{3}$$

The correct way to do this is write
$$y[n] = (x\star h)[n] = \sum_{k=-\infty}^\infty x[k]h[n-k]$$
and then substitute $n-k$ for $n$ in $h[n]=u[n]-u[n-8]$
to
get $h[n-k] = u[n-k]- u[n-k-8]$ and not $u[-n-k]-u[(-n-8)-k]$
the way that the OP has it.

Some initial thought might have
shown that (assuming that $u[n]$ is the unit step),
that $h[n]$ is $0$ for $n < 0$ (since then both $u[n]$ and $u[n-8]$
are $0$, and that $h[n] = 0$ for $n \geq 8$ (since then
both $u[n]$ and $u[n-8]$ are $1$. Consequently for any given
fixed integer $n$, $h[n-k]$ is nonzero only for
$$0 \leq n-k \leq 7 \Rightarrow n-7 \leq k \leq n$$
and so the limits on the sum could be reduced from
$-\infty$ and $\infty$ to $n-7$ and $n$ (cf. my comment
on the OP's question).  I will leave it to the OP to figure
out whether the fact that $x[n]$ is also time-limited
can be used to further reduce the limits on the sum.
Hint: the answer to this might depend on the specific
value of $n$, and you may need to try several different
values of $n$, e.g. $0, 3, 8, 10, -40$ etc to see if you
can discern a pattern.


applicability of dictionary learning in compressed sensing problem Compressed Sensing Problem:

$Y = MX$, $M$ = measurement matrix (known), $X$ = full signal (unknown), $Y$ = sampled points (known). Objective is to obtain $X$ using the concept of sparsity i.e. $X = \phi S$, $S$ = sparse vector (unknown), and $\phi$ = a basis which converts the signal to a sparse domain (known). So $Y=M \phi S$, from this experssion we find $S$, and put it into $X = \phi S$ to get $X$.

Dictionary Learning Problem

$P = \psi Q$, $P$ is a signal (known), $Q$ = sparse vector (unknown), and $\psi$ = a basis which converts the signal $P$ to a sparse domain (unknown). We find $\psi$ and $Q$ using alternating minimization algorithms like k-svd.

My question is: Compressed sensing using Dictionary Learning

$Y=M \phi S$, $Y$ = sampled points (known), $M$ = measurement matrix (known), $S$ = sparse vector (unknown), and $\phi$ = a basis which converts the signal $Y$ to a sparse domain (unknown). Is this possible to solve at all?

Referring to appropriate references will also help a lot.

Why of course, it is certainly possible. As you already noted, dictionary learning is tightly coupled to sparse recovery, which is what you need for compressed sensing.

In DL, we are trying to factor $Y = D X$, which is only possible if additional assumptions are made. Typically that includes some notion of sparsity of $X$ and some additional normalization constraints on $D$. There are a few more conditions though to even render the problem unique. Note that $D$ might not be a basis, in general it is not even square. Typical strategies for DL include alternating optimization (have a look for MOD-based algorithms): starting from some guess about one of them, estimate the other one, and then take turns. This takes advantage of the fact that while finding both $D$ and $X$ is a hard (non-convex) problem, finding one of them given the other is much easier:

Finding $D$ given $X$ (and $Y$) is a linear least squares problem
Finding $X$ given $D$ (and $Y$) is a sparse recovery problem (which by itself is also non-convex but has nice convex relaxations that are quite well studied).

In the CS setup, we are encountering a similar situation, except that here, $D$ is typically composed of a measurement matrix we know (say, $\Phi$) and some sparsifying dictionary or basis (say, $A$) we may or may not know so that $D=\Phi A$ and $Y=D X = \Phi A X$. In this case $A$ could be a basis (i.e., it is square) but it could also be some more general dictionary (i.e., flat). The matrix $\Phi$ on the other hand is always flat. Given that, the above strategy still applies and you can apply alternating optimization:

Finding $A$ given $\Phi$, $X$, and $Y$ is a linear least squares problem (unless you need $A$ to be a basis, which adds an orthonormality constraint)
Finding $X$ given $\Phi$, $A$, and $Y$ is a sparse recovery problem (in no way different from what is already used in the CS context).

That said, there are strategies to optimize $D$ and $X$ jointly, like the K-SVD. They might not be directly applicable.


How to find Fourier Series Coefficients I saw many solved examples about this topic but again I coudn't come up with any solutions about this question. How can I find the Fourier Series coefficients of the following signal ?

$x(t)=2 \cos(3\pi t) + \sin(100\pi t+\frac{\pi }{3})$

I know that; $\cos(θ)=\frac12 (e^{jθ} + e^{−jθ})$  and,
$\sin(θ)=\frac{1}{2j} (e^{jθ} - e^{−jθ})$

I also know that I should use
$$x(t) = \sum_{-\infty}^{\infty} a_{k} e^{jk(2\pi/T)t}$$

But I'm having trouble to define the fundamental period $T$ and the relation between sinusoidal terms and coefficients $a_k$, to sum all things together. Thanks to everyone...

First, recognise that the fundamental period of the signal x(t) (sum of two sinusoidals with periods 2/3 s and 1/50 s) is 2s. And the fundamental frequency, $f_0$, is 0.5 hz.

Then considering the CT Fourier Series expansion of the periodic signal x(t) wrt its fundamental frequency $f_0$ given as $\sum_{-\infty}^{\infty}{a_k e^{jk2\pi f_0t}}$ you can  therefore find the coefficients corresponding to $2\cos(3\pi t)$ and $\sin(100\pi t + \pi/3)$ with k=3 and k=100 respectively.

Then, by expressing the sinusoidals in terms of complex exponentials using the Euler's identity this becomes $a_3=1 ,~ a_{-3}= 1 , ~~a_{100} = e^{j\pi/3}/2j, ~~a_{-100} = -e^{-j\pi/3}/2j $  which can be further manipulated if you want a rectengular representation of the complex coefficient instead of the above mixed-polar one.

Assuming no errors made during the arithmetic.


Why do these 2 methods give different solutions? I need to solve what is underlined in red for $x_i$, nut currently I'm interested in the right side of the equation only.

On the left I sarted by doing the Laplace transform of $x_u'$ and $x_u$, and this gives the correct solution.

On the right side I tried to derive $x_u$ first, plug everything into the equation and then do the Laplace transform, and this gives a completely different solution. Is this method incorrect? Or am I doing something wrong?

The problem is that you took the derivative of the function

$$\hat{x}_u(t)=2e^{-3t}-e^{-4t}\tag{1}$$

whereas using the Laplace transform you implicitly assumed that $x_u(t)$ equals zero for $t<0$:

$$x_u(t)=\hat{x}_u(t)u(t)=(2e^{-3t}-e^{-4t})u(t)\tag{2}$$

where $u(t)$ is the unit step function.

If you take the derivative of $(2)$ then you'll get the same result as with the Laplace transform, taking into account that the derivative of the step function $u(t)$ is the Dirac delta impulse $\delta(t)$:

$$\frac{dx_u(t)}{dt}=\frac{d\hat{x}_u(t)}{dt}u(t)+\hat{x}_u(t)\delta(t)=\frac{d\hat{x}_u(t)}{dt}u(t)+\hat{x}_u(0)\delta(t)\tag{3}$$

where I've used the fact that for any function $f(t)$ that is continuous at $t=0$ we have $f(t)\delta(t)=f(0)\delta(t)$. I trust that you can take it from here.


Complex exp. Fourier series, finding $x(t)$ when $X(j\omega)$ is given as magnitude and phase plot I'm watching Neso Academy series on Signals and Systems, and in one of the videos the problem is to find $x(t)$ when magnitude and phase plot are given. The plot looks like this:

When he finishes calculation he gets:

$4 + 4cos(3 \omega t+ \dfrac\pi2) + 8cos(4 \omega t- \dfrac\pi2)$

I understand the steps that are required to get this, but I don't understand why is it cosine? When you look at magnitude plot you have just 3 components: DC offset = 4 and 2 harmonics.

As far as I know the values on magnitude plot represent sine value and phase plot represents just phase shift of that sine wave. Without doing anything but looking at the plot I would write following:

$4 + 4sin (3 \omega t+ \dfrac\pi2) + 8 sin(4wt-\dfrac\pi2)$

And it's completely wrong, why is it so?

BTW link to the video.

The video could've stopped at the answer:

$x(t)=4+ 2e^{j(3\omega t + \pi /2)} + 4e^{j(4\omega t - \pi /2)} + 2e^{-j(3\omega t + \pi /2)} + 4e^{-j(4\omega t - \pi /2)}$

And this answer can be read directly from the plots. The idea is that the $n^{th}$ term in the sum is equal to $|C_n|e^{j(n\omega t+\angle C_n)}$. You get this from the Fourier series equation which is $x(t)=\sum_n C_n e^{jn\omega t}$, where $C_n$ may be a complex coefficient meaning it has a magnitude and a phase. Re-write $C_n$ as $C_n = |C_n|e^{j\angle C_n}$, plug into the Fourier series equation and you get what I wrote above.

The extra steps that the video does is just to convert the complex exponential expression into cos/sin functions using Euler's formula. One result of Euler's formula is that $2cos(\theta)=e^{j\theta}+e^{-j\theta}$, and that is why in the end it turns out to be cosine and not sine.


Can we see the Fourier transform as a filtering operation? Is the Fourier Transform (FT) a filtering operation?
Is it possible to graphically represent this transform?
I know that a signal can be represented in the frequency domain, but I want to know if the FT operation can be "drawn" as a coefficient series.

The FT matrix is composed by $N\times N$ elements, and I can apply it to the input sequence of length $N$ in order to obtain an output sequence of length $N$. But I want to understand if there is another way to see the FT matrix in one dimension.

Hint:

The DFT matrix only contains powers of a complex root of unit. Every row or column holds increasing powers of a common power ($\omega^{mn}=(\omega^m)^n$); hence the real and imaginary parts follow (co)sinusoids of increasing frequencies.

https://www.dsprelated.com/freebooks/mdft/DFT_Sinusoids.html

Every component of the DFT of a signal is a convolution of the signal with those sinusoids.


Plotting magnitude and phase for frequency spectrum I would like to find the frequency spectrum by drawing magnitude and phase for the following signal:

S=sin(2*pi*100*t)+cos(pi*500*t);

X=sinc(2*pi*t);

Is there something particular to this problem that I should be noticing?

We can easily find the spectra giving a time domain like

t=0:0.01:0.1;

And using the fft command for both signals

X=fft(S);

What did they mean by drawing the magnitude and phase?

for plotting magnitude spectrum versus frequency you can use the following script:

t = 0:0.001:0.2;
s = sin(2*pi*100*t)+cos(pi*500*t);
S = abs(fftshift(fft(s)));

Fs = 1000;
df = 1000/(length(s));
f = -Fs/2:df:Fs/2-df;

plot(f,S);

Here I used your first signal and plotted its magnitude spectrum. You must notice for your signal min choosable Nyquist frequency is 500 Hz But we use it two times of it(which is 1000Hz), because we also want to see the part before f =0 in the frequency spectrum. Using this we find df(delta f) frequency spacing. frequency vector must be multiples of df. By this way spectrum will be drawn correctly.

Here is the output:


Why is each window/frame overlapping? I am currently reading up on some speech recognition theory.
As it speech signals usually are non stationary, but if a speech is seen in a very small time frame can it be deemed stationary, as would only consist of the sound of phoneme.

As explained: http://recognize-speech.com/preprocessing#[object HTMLHeadingElement]

But then this comes up:

  As it can be seen in the figure above there can be an overlapping of
  the windows. The frame length and window length is dependent on the
  scope of application and algorithms used. In speech processing the
  value for the frame length  typically varies between 5 to 25ms and for
  the window length  between 20 and 25ms [ref.]. Smaller overlapping
  means larger time shift in the signal, therefore lower processor
  demand, but the difference of parameter values (e.g. feature vectors)
  of neighbouring frames can be higher. Whereas larger overlapping can
  result in a smoother change of the parameter values of the frames,
  although higher processing power is needed.

My question is how the overlap lower/increase the processor demand?

Why is each window/frame overlapping?

Windowing is a means to stationarize signals. Inside a small enough window, you can expect that the  properties of the signal chunk do not vary too fast. And you now can use tools well-suited to stationary signals, like Fourier-based techniques.
You can imagine non-overlapping rectangular windows, each defining a frame. Each sample is, somehow, treated with the same weight (the height of the window). However, when you look at the features extracted from two consecutive frames, or when processing them, the change of property between the frames may induce a discontinuity, or a jump ("the difference of parameter values  of neighbouring frames can be higher"). The same phenomenon happens with JPEG images, using $8\times 8$ non-overlapping blocks. This tiling or blocking effect is already annoying for the eye, it is even more disturbing for the ear.

Adding to this, if the signal evolves slowly with time, the left-end or the right-end sample could be seen as the "most different" samples: the frame is inherent localized at a specific location, depending on the window shape. If the window is symmetric, the central location is the center of the frame (some use asymmetric windows for speech processing).

So you can imagine giving some smaller weights to the edges of the frame. This is what you do when you are averaging a signal with unimodal weights, like triangular or Gaussian, instead of uniform. This can be beneficial for frequency estimation, as other windows have better behavior than the uniform one. However, this creates some imbalance in the processing (even information loss as said by @hotpaw2, esp. if the window vanishes) of some samples (the edge ones), unless you compensate it in some way, for instance by choosing a sum of windows whose amplitude ("constant overlap-add") or energy is close to one. Some figures can be found in
Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new at-top windows, and the framework of oversampled filter banks provides a quite unified treatment of such constraints.

My question is how the overlap lower/increase the processor demand?

It depends a lot of the implementation.  When the processing is done offline, you can save a lot of redundant operations, by folding the windows into independent blocks, subsampling where you can, etc. Using polyphase matrices or lifting from the overlapped filter bank theory can help you a lot in saving operations. Also, CPU and GPU aware processing can completely change the computational burden deal.
So what follows is quite crude.

Without much optimization however, if you have an overlap of $p$%, each sample will be processed, in average, $c = 1/(1-p)$ times. This is about the factor for the  number of frames that could be naively processed: $c=2$ for $p=50\%$, $c=4$ for $p=75\%$. If $p= 0$, frames can easily be processed in parallel. Otherwise, some memory or processing is required for the overlap.

At the extreme, you can perform a windowing centered at each sample. This is the principle at the basis of the time-frequency analysis. The computational cost can be high. But the high redundancy has a good side: any processing in the frequency domain would be more robust to noise, or imperfections in the processing: you should hear less clicks and pops.


Relation between a Fourier series harmonic component and its power My question is about the meaning of power spectrum derived from the Fourier series coefficients.

Fourier series is shown below:

$$f(t)=a_0+\sum_{n=1}^{\infty}a_n\sin(\omega_n t)+\sum_{n=1}^{\infty}b_n\cos(\omega_n t) $$

where
$\omega_n = 2\pi n f $

Since the Fourier-series above is summing sines and cosines for all the harmonic frequencies to obtain a periodic signal $f(t)$, at a particular frequency the the sine and cosine sum together with offset can be written as:

$$a_0 + a_n\sin(\omega_n t) + b_n \cos(\omega_n t) $$

One can show the Fourier-series coefficient for each harmonic as:

But textbooks shows the power spectrum instead as:

Why is the power magnitude equal to the square root of the sum of squares of the sine and cosine components? How can I derive this:

$$\text{Magnitude} = \sqrt{a_n^2 + b_n^2}$$

from the following:

$$ a_0 + a_n\sin(\omega_n t) + b_n \cos(\omega_n t) $$

See first of all, you have confused magnitude spectrum to power spectrum. The figure you have shown in green is not power but amplitude spectrum hence y axis is represented as "Magnitude"

what is happening in the figure you have shown is the periodic signal is expressed as a trigonometric series where

$$x(t)= a_0 +\sum_{k=1}^{\infty} [a_k\cos(k\omega_0t) +b_k\sin(k\omega_0t)]$$

where $a_k$ and $b_k$ are trigonometric Fourier series coefficients and $\omega_0$ is fundamental period of signal x(t).

Now to express the amplitude and phase of the periodic signal in terms of Fourier series coefficients, trigonometric series representation is inconvenient as determining phase and magnitude of spectrum cannot be done directly from the above mentioned equation.

So we generally express periodic signal in terms of polar form of Fourier series where the polar form is given by
$$x(t) = M_0 + \sum_{k=1}^{\infty} M_k \cos(k\omega_0t - \phi_k)$$

where from the above expression you can see that $M_k$ denotes amplitude at kth harmonic and $\phi_k$ represents the phase of the kth harmonic.

where
$$M_0 =a_0$$
To find the remaining terms $a_k$ and $b_k$ in terms of $M_k$ and $\phi_k$, let's expand the summation of both the expressions of $x(t)$ and equate them.

i.e.  $$a_k\cos(k\omega_0t) +b_k\sin(k\omega_0t) = M_k\cos(k\omega_0t-\phi_k)$$

now expand the right hand side equation which is of form $\cos(A-B)$

So $$a_k \cos(k\omega_0t)+b_k \sin(k\omega_0t) = M_k\cos(k\omega_0t)\cos(\phi_k) +M_k \sin(k\omega_0t)\sin(\phi_k)$$

so equating coefficients on both sides we get

$$ a_k = M_k \cos(\phi_k)$$

and $$b_k = M_k \sin(\phi_k)$$

now to find magnitude $M_k$, square and add both the above terms, then we get

$$M_k =\sqrt{a_k^2 +b_k^2}$$

and phase is obtained by dividing the expressions for $a_k$ and $b_k$

$$\phi_k = \arctan(b_k/a_k)$$

Now the graph you have shown in green is the plotted $M_k$ values.

Note that you can also find magnitude and phase using exponential Fourier series coefficients, but then you will have both positive and negative harmonics in exponential Fourier series unlike trigonometric Fourier series which has only positive harmonics. So polar form can be seen as a direct equivalence with the trigonometric Fourier series coefficients.


which programming language is the best for image processing? I used matlab for my projects but wondering is there any framework similar but easier than matlab.

Most image processing code in image signal processors
(ISP) in the industry are written in C++ and assembly, which I would say have more syntax and language gibberish than MATLAB.

So if you are considering to write code for educational purposes than MATLAB and python have good libraries and easier to code.

Else if you are writing code to be run on real HW, C++ and assembly mostly.

Open CV is a good cross platform library with C++ and other languages supported as well


What's a good physical description of the following matlab code? I've inherited some matlab code and I can't really understand what the resultant signal really means. The input is audio and the code first performs an STFT/spectrogram (sliding window FFT of size 1024 with 50% overlap, hanning window applied). The output is then filtered to the range (using Matlab notation) to the range (2:N/2). In other words, the first bin is dropped and the negative half frequencies too. Just to be clear, for each 1024 window of audio data (audio_chunk below) the code does:

y = abs(fft(hann_win.*audio_chunk));
filt_y = y(2:length(y)/2);

The value filt_y is then summed and that becomes the output signal:

output = sum(filt_y);

So, for every window's FFT we get a single output value. The output sampling frequency is audio_sampling_frequency/512 because we slide the 1024 sample window along by 512 samples at each step.

My question is: what would you call that output signal? What sort of physical significance does it have, or what use is it for analysis?

The input is audio and the code first performs an STFT/spectrogram (sliding window FFT of size 1024 with 50% overlap, hanning window applied). The output is then filtered to the range (using Matlab notation) to the range (2:N/2). In other words, the first bin is dropped and the negative half frequencies too.
. . .
The value filt_y is then summed and that becomes the output signal:

What you are describing thus far is equivalent to a boxcar bandpass filter followed by integration similar to the way a moving average filter works.

So, for every window's FFT we get a single output value. The output frequency is audio_frequency/512 because we slide the 1024 sample window along by 512 samples at each step.

It would not be straightforward to define "output frequency" here. Most of the times, when you refer to "frequency" you might implicitly be referring to sinusoids. However, here, there is no sinusoidal output. But, because of properties of the DFT, your output will indeed fluctuate at a frequency that depends on the length of the FFT window.
This is because, the FFT assumes that the signal at its input is periodic. In other words, it assumes that whatever lands in its input frame is played again and again in a loop to infinity on either sides (both positive and negative). But this is not what is happening in reality. Because of this, your "running sums" from frame to frame have discontinuities. If you were to listen to that signal, IT!WOU!'D.CO!NTAI!NDIS!CON!TIN!UIT!IES!.
There are however two elements in your code that reduce the effect of those discontinuities. The first, is the use of the hann window and the second is the fact that the summation drops the DC component (and low frequencies depending on NFFT). With the use of the window, the discontinuities between the end of the frame and its beginning are reduced (remember, that is the start and end of the "loop"). By dropping the DC, you drop the result of partial summation that could end up having very high values.
Imagine for example that you are summing a zero mean sinusoid with a "frame" whose length just so happens to be half its period. Then in that case, during the positive phase of the sinusoid you would get a high positive value and then during the negative phase you would get a negative value. If the length of the "frame" was one period of the sinusoid, you would get a zero mean.
If you left the DC in and tried to "listen" to the signal you would be getting a soft "flutter" sound in the background, like holding a microphone close to a fan because of these discontinuities that are produced at each frame.

My question is: what would you call that output signal?

Some sort of band-limited integration. As mentioned previously, this is an integrator with a crude band-pass filter at its input.

What sort of physical significance does it have, or what use is it for analysis?

It returns the total amplitude of the frame of the signal minus its average value. It would be related to the average amplitude in a "moving average filter" sense.
It is difficult to say what its use would be for analysis without knowing what else is going on in the processing pipeline, but something like this could be used to get a sense of how "loud" is the current frame.
If this was a compressor for example, this signal could be used to trigger it at a sharp transition from soft sound to louder sound. In fact, because of the crude band-pass filter it contains at its input, this would then become a de-esser.
Hope this helps.


Do lowpass filters affect the integral over the signal? For which filters, does the resulting signal have the same integral (/average) as the original? Obviously this is true for a moving average, but is it also mathematically true for e.g. a Butterworth filter?

I am trying to build a lowpass filter for a signal of measurement samples. While there is a lot of noise, the average (over a sizable number of samples) is correct. I need to filter it retaining that property, regardless of the input frequency spectrum.

You need a filter whose transfer function evaluated at 0 (DC component) is 1, or equivalently, whose impulse response sums to 1.


hough voting, turning transformation into a vote I need help understanding how/what does it mean to "turn" a transformation parameters into a vote.
I know how to recover the transformation parameters (in my case: 6 DOF), but how this turns into a vote in 6D (??) space ?

Thank You !

Determine the accuracy you will need for each of your six DOF. Don't use larger accuracy than you actually need, since it will get you in trouble later, when you have to find maxima in the Hough Space.
Determine the range of possible transformation parameters for you six DOF.
Now create a six dimensional array A (one dimension for each DOF) where each dimension has the size range/accuracy for the current DOF. This array is your Hough Space.
Now, for any given vector v = [DOF1, DOF2, DOF3, DOF4, DOF5, DOF6] increment the bin in A in which v falls. You have now turned this single vector of transformation parameters onto a vote in the Hough Space,

2D example:

x is in the range [1;6] and y also in the range [1;6].
x should be found with an accuracy of one whereas y only needs an accuracy of two.
A(x,y) looks like this:

0;0;0
0;0;0
0;0;0
0;0;0
0;0;0
0;0;0

where the bold bin corresponds to x=[3] and y=[1;2].

So if you would like to cast a vote for a vector v(x,y)=[3,1] you would simply increment the bold zero by one. Votes for v(x,y)=[3,2] would in this case fall in the same bin, because of the reduced accuracy needed for y, whereas votes for v(x,y)=[3,3] would go into the bin to the right.

I hope this was helpful, and did not just confuse you more. ;)


Efficient correlation of a low duty cycle training sequence Is there a way to efficiently correlate a training sequence that is N samples long, framed at M samples where M >> N, with L occurrences of such frames (see below). For the pedants, the training sequence may or may not vary in the L subsequent frames, thus the longer correlation template.

The samples in between the subsequent training sequences would be unknown so they would be 0's in the correlation template. There is a lot of wasted processing with the multiplications/additions of all the 0's. I am wondering what optimizations can be made by taking advantage of such a signal structure.

Example:

       <- N ->
        _____
       |     |         ... repeated L-1 times ...
 ______|     |_______
 <------- M --------->

One way to approach this would be to have a bank of $L$ subcorrelators, each of which use the corresponding $N$-sample long template. Then, in order to get the effect of correlating against the longer gated sequence, one need simply to sum the outputs of the subcorrelators with the appropriate delays applied. This is equivalent to breaking the long matched filter into a parallel formulation with $L$ subfilters.

For instance, assume the output of subcorrelator $i$ is $c_i[n]$, and that the notional output of the full $ML$-sample corelator is $c[n]$. The full correlator output is just equal to the sum of the shorter subcorrelators, after taking into account the $M$-sample delay between the expected position of each template in the input signal:

$$
c[n] = \sum_{i=0}^{L-1} c_i[n-(L-1-i)M]
$$

Using this structure, for each output sample, you effectively implement an $LN$-tap FIR filter instead of an $LM$-tap one. With that said, if any of these numbers are large (greater than 50-60 or so), it's likely more efficient to use FFT-based fast convolution techniques instead. This method is compatible with that; you would simply use fast convolution to calculate each of the short subcorrelator outputs, then perform the delayed sum as described above. I've implemented this very technique before with good success.

There are a couple other things to think about with this approach:

The OP suggests that some of the $L$ frames may not have distinct patterns from each other. If this is the case, you can achieve some complexity reduction by only implementing the number of unique subpatterns, then tapping the outputs from each subcorrelator appropriately.
This sort of technique can run into issues with scalloping loss if the overall length of the correlator ($LM$ samples) is large with respect to any carrier frequency uncertainty on the input signal. This is inherent with any kind of coherent correlator structure; the longer you correlate (thus giving you more potential gain), the more frequency selective the correlator is (the more loss you have for a given frequency offset in your signal).


image compression versus image denoising Wavelet based image denoising may be performed by thresholding (selecting a threshold value, and discarding all values below the threshold. Wavelet-based image compression  may also be performed by a similar way. What is the difference then between the two operations? can we consider denoising a certain type of compression?
Discarding in case of denoising means removing the coefficients below the threshold, while in compression all values below the threshold are set to zero? why setting the values to zero reduces the entropy? they are not removed, and if the image is transmitted across a wireless channel, the zero values are coded and transmitted.

Perceptually lossless (but actually lossy) image compression aims to maintain visually the same image while reducing the number of bits necessary to represent it. If you try to compress aggressively however, visually apparent compression artifacts will show up.

Wavelet based image denoising aims to change the appearance of an image by getting rid of what is considered as noise. We expect the denoised image to look better, compared to the noisy version.

One side effect of image denoising is also a reduction in the size of the image, as you are effectievely discarding some data from the image. But it will never be as efficient as a compression algorithm.

Furthermore, since noise is always considered as irrelevant information, it can be useful to denoise an image (if necessary) before compressing it, allowing you to achieve higher compression rations.

However, note that image denoisers are not perfect algorithms in terms of discriminating noise from true image data, and they may and will discard not only noise but also some very critical image information (especially at high frequencies, textures an edges).


What the disadvantages of up-sampling Let's say we have a signal $x[n]$ upsampled at the transmitter by $4$ and then downsampled at the receiver by $4$ too.   As known, the performance will be much better compared with transmitting the signal without upsampling.  but what I need to understand, what is the disadvantage of using up-sampling ? Does that effect the data rate of the transmitted signal since we transmit $4$ times the signal $x[n]$ instead of $x[n]$?

When you upsample a signal by a factor of 4, then, you have to churn out 4 times the number of samples in the same time in order to maintain the same Symbol Rate. Otherwise your symbol rate will be $\frac{1}{4}$ times the original symbol rate. This means your DAC needs to operate at least at 4 times the sampling frequency now. This can be costly. Also faster clocks need better synchronization techniques.


Investigating correlation for unequal signal lengths I have measurements from two tablet sensors. The first sensor measures the touch pressure applied to the tablet and the second sensor is the acceleration of the tablet (in z direction). The sensor data are stored in csv files where the first column contains the timestamp (in milliseconds) and the second column contains the sensor data (either touch pressure or acceleration). I have uploaded the two sensors files:
https://1drv.ms/f/s!AlkQuYgB1McYgudFUrpMzMopUr2wVw

Now I would like to investigate if there exists a correlation between touch pressure and acceleration. I could just calculate the cross-correlation or plotting the pressure against the acceleration (and perhaps add a linear regression line).

The problem is that the acceleration measurements have a more or less constant sampling rate but for the pressure I have not a constant sampling rate. That means that the number of measurements for acceleration and pressure is also not the same (and the timestamps are therefore also not the same).

How can I still calculate the correlation and producing a accleration vs. pressure plot? I thought about a interpolation of the accleration measurements to the timestamps of the pressure readings but I don't know what type of interpolation to use (linear or cubic seems not to be right). Perhaps there is a better way of doing it.

Here are the plots of the two sensors:

If you are using MATLAB, there is a built in object type, called

  timeseries

, that is very useful for this sort of problem.

>    Methods for class timeseries:
>
> addevent              getabstime            gettsbeforeatevent    mean
> pvset                 times                  addsample
> getdatasamples        gettsbeforeevent      median
> rdivide               timeseries             append
> getdatasamplesize     gettsbetweenevents    min
> resample              tsprops                createTstoolNode
> getinterpmethod       hasduplicatetimes     minus                 set
> utArithCommonData      delevent              getprop
> idealfilter           mldivide              setabstime
> utArithCommonOutput    delsample             getqualitydesc
> init                  mode                  setinterpmethod
> utArithCommonTime      detrend               getsamples            iqr
> mrdivide              setprop               utChkforSlashInName
> display               getsampleusingtime    isequal
> mtimes                setuniformtime        utGetEventTime         eq
> gettsafteratevent     isequalwithequalnans  plot                  std
> utStatCalculation      filter                gettsafterevent
> ldivide               plus                  sum                   var
> get                   gettsatevent          max
> pvget                 synchronize
>
> Static methods:
>
> loadobj               tsChkTime             tsgetrelativetime

The setuniformtime method will do the interpolation for you. The other methods are often useful.