How to represent impulse function in 2D? To be more specific I want to show that impulse function in 2D can be represented as $β(r)=δ(r)/πr$.

Also I want to show that each projection of a two dimensional impulse function at the origin is a Delta function.

I could not find any useful information about these two problems, if anyone could help I will be appreciated.

So dealing with generalized functions like the Dirac delta requires some care, and when dealing with N-dimensional versions you need to be very explicit with your notation to keep things straight.

I'll denote the 2 dimensional delta function in polar coordinates at the origin as ${}^2\delta(r, \theta) = {}^2\delta(r)$, since for the special case of the origin, $\theta$ doesn't matter.

For this derivation, ${}^2\delta(r)$ represents the following limiting sequence of functions (that are asymmetrical about 0):

$${}^2\delta(r) = \lim_{\epsilon \rightarrow 0^+} \dfrac{1}{\pi\epsilon^2} \quad 0 < r < \epsilon$$

Using this limiting sequence of functions, the 2-D delta function at the origin in polar coordinates has the following property of "integrating to 1" for the 2-D integration:

$$\int_0^{2\pi} \int_0^{\infty} {}^2\delta(r) \space r \space \mathrm{d}r \space \mathrm{d}\theta = 1$$

To separate the $r$ and $\theta$ portions of ${}^2\delta(r)$, so we can use a 1D Dirac Delta function, I'll use the following limiting sequence of functions (that are asymmetric about 0) for the 1D Dirac Delta function:

$${}^1\delta_a(r) = \lim_{\epsilon \rightarrow 0^+} \dfrac{1}{\epsilon} \quad 0 < r < \epsilon$$

Using this limiting sequence of functions, the 1-D delta function at the origin has the following property of "integrating to 1" for the 1-D integration:
$$\int_0^{\infty} {}^1\delta_a(r) \space \mathrm{d}r = 1$$

We can then equate the two above integrals to get the relationship between  $ {}^2\delta(r)$ and ${}^1\delta_a(r)$

$$\begin{align}\\
\int_0^{\infty} {}^1\delta_a(r) \space \mathrm{d}r &= 1\\
&= \int_0^{2\pi} \int_0^{\infty} {}^2\delta(r) \space r \space \mathrm{d}r \space \mathrm{d}\theta\\
&= \int_0^{\infty} {}^2\delta(r) \space 2\pi r \space \mathrm{d}r\\
\end{align}$$

so we have

$${}^2\delta(r) = \dfrac{{}^1\delta_a(r)}{2\pi r}$$

This differs from your desired result by a factor of $2$, because I did not allow the 1 dimensional delta function to be symmertic around $0$, which wouldn't make sense for a polar origin where $r$ cannot be less than $0$ and $\theta$ ranges in $[0, 2\pi)$.

If one allows $r$ to go negative and $\theta$ to range in $[0, \pi)$, then one can show for a 1 dimensional delta function that is symmetric around $0$
$${}^2\delta(r) = \dfrac{{}^1\delta(r)}{\pi |r|}$$

(note the absolute value bars!)

This is a good example of why delta functions require careful definitions and considerations when making statements.

To show the projections of a 2-D delta function are themselves delta functions, start with a suitable limiting sequence of functions for ${}^2\delta(x,y)$, like an infinitely thin and tall rectangle function both in x and in y, such that

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} {}^2\delta(x,y) \space \mathrm{d}x \space \mathrm{d}y = 1$$

and the answer should fall out from the limiting sequence of functions you use.


Random sampling vs uniform sampling In this paper of Lustig, he speaks about a something which appears unintuitive: sampling at random may exhibit better performance than sampling uniformly. I tried to understand this starting from page 15 of these slides, but I can't really make sense of anything.

Why, if we take random permutation of frequency coefficients, do we get a better reconstruction in terms of signal similarity? Why does this give better reconstruction, and what's the intuition behind the phenomenon?

The key idea is that the random sampling approach enforces more constraints on the resulting signal than the uniform sampling approach does.

The POCS (projections onto convex sets) algorithm used for the reconstruction of the randomly sampled signal is the key piece: it enforces:

That the signal must be from this spectrum.
That the signal is real-valued.
What we know about the signal's spectrum (the randomly sampled Fourier coefficients).
What we know about the signal's time-domain form.

The uniform sampling approach makes no attempt to enforce the bold constraint.

Here's an example which shows:

Top left: Uniform undersampling and reconstruction using the FFT only.
Top right: Random undersampling and reconstruction using the FFT only.
Bottom left: Reconstruction using POCS

As you can see, that final constrain really improves the reconstruction hugely.

R Code Below

#30219

T <- 128
N <- 5

x <- rep(0, T)
x[sample(T,N)] <- rep(1,N)

X <- fft(x);
Xu <- rep(0, T)
Xu[seq(1,T,4)] <- X[seq(1,T,4)];
xu <- fft(Xu, inverse = TRUE)*4/T;

par(mfrow=c(2,2))
plot(x, type="l", lwd = 5, ylim = c(0,1.2))
lines(abs(xu), col="red")
title('Original (black) & reconstructed\n from uniform undersampling (red)')

Xr <- rep(0,T)
r_ix <- sample(T,T/4)
Xr[r_ix] <- X[r_ix]
xr <- fft(Xr, inverse = TRUE)*4/T

plot(x, type="l", lwd = 5, ylim = c(0,1.2))
lines(abs(xr), col="red")
#lines(Re(xr), col="blue")
#lines(Im(xr), col="green")
title('Original (black) & reconstructed\n from non-uniform undersampling (red)')

#soft thresh function

softThresh <- function(vals_to_threshold, lambda)
{
  ix <- which(abs(vals_to_threshold) < lambda)
  vals_to_threshold[ix] <- rep(0, length(ix))

  ix <- which(vals_to_threshold >= lambda)
  vals_to_threshold[ix] <- vals_to_threshold[ix] - lambda

  ix <- which(vals_to_threshold <= -lambda)
  vals_to_threshold[ix] <- vals_to_threshold[ix] + lambda

  return(vals_to_threshold)
}

# POCS
lambda <- 0.1
Xhat <- Xr
for (iteration in seq(1,100))
{
  # 1. Compute the inverse FT to get estimate
  xhat <- Re(fft(Xhat, inverse = TRUE)/T)
  # 2. Apply Softrhesh in the time domain
  xhat <- softThresh(xhat, lambda)
  # 3. Find the FFT
  Xhat <- fft(xhat)
  # 4. Enforce known values
  Xhat[r_ix] <- X[r_ix]
}

plot(x, type="l", lwd = 5, ylim = c(0,1.2))
lines(xhat, col="red")
title('Reconstructed using POCS')


Missing harmonics (power-line) I have some biomedical signals. If I calculate their magnitude spectrums, I see power-line interference at 50 Hz. I was excepting a harmonic at 100 Hz, but it turns out there is none. However, there is a harmonic at 150 Hz. Is this common or should I worry that I have done something wrong?

Yes, this is common. Common non-linearities (such as op-amp clipping) preserve the symmetry of the upper and lower section of the waveform; and as such, they do not introduce even harmonics - just odd harmonics.


DFT for audio classification - signals of different lengths If I have a set of signals each with $n$ points, I can take the DFT of each and then use the resulting frequency domain vectors as features for some classification algorithm. If the signals are differing lengths, the frequency basis used by the DFT is no longer the same. This means I can no longer directly use DFT as features if I want to compare signals of varying length. What should I do about this?

In other words, I want to convert each audio signal into a length $m$ vector of frequency components, and use the same $m$ for all signals.

I realise there is fundamentally more information in a longer signal, so any conversion will involve some loss of information.

As you explicitly mention, that you are aware of a loss of information and are willing to accept that, a possible solution will be modulo-N reduction or wrapping, as described by Orfanidis in the freely available book "Introduction to Signal Processing" on pages 489ff.

In modulo-$N$ reduction, a vector $x$ with length $L > N$, where $N$ is the DFT length, is divided into non-overlapping blocks of length $N$. These blocks are summed up, resulting in one vector $\tilde{x}$ of length $N$ .

Figure taken from Orfanidis (linked above), Fig. 9.5.1 on page 489.

You then apply the DFT on $\tilde{x}$, which results in $\tilde{X} = \text{DFT}\{ \tilde{x} \}$ with length $N$. Orfanidis proves that $\tilde{X}$ is exactly equal to the $N$-point DFT of $x$, so is equal to $X$.

In your case, that would mean you could find a minimal length $N$, which is guaranteed for all input signals. You then apply a modulo-$N$ reduction on all your input signals and compare these length-$N$ DFTs. Of course, a higher resolution would be possible for longer signals, and you just ignore this by using the modulo-$N$ reduction method.

Similarly, you could just zero-pad all signals to a maximal length $M$. All your signals will then have a high frequency resolution defined by $M$. Important point: you do have that computational frequency resolution, but the physical frequency resolution is still determined by the length of your signals $L$. More details on this are given in any decent signal processing class or on pages 482ff in Orfanidis.


Difference Between Correlation and Convolution in the Context of Image Processing Could you please clearly explain what is the difference between correlation and convolution that is done by a filter on an image?

I mean in terms of signal processing definition I know that convolution describes the output of an LTI system, that is if an LTI system produces an output due to convolution with an input system then the output signal can be described as the result of convolution of the input signal and the impulse response of the LTI system. As for the correlation, it describes the similarities between to signals. But how does convolution and correlation effect on a image and how different are they in terms of effects?

Thanks

Convolution is correlation with the filter rotated 180 degrees. This makes no difference, if the filter is symmetric, like a Gaussian, or a Laplacian.  But it makes a whole lot of difference, when the filter is not symmetric, like a derivative.

The reason we need convolution is that it is associative, while correlation, in general, is not. To see why this is true, remember that convolution is multiplication in the frequency domain, which is obviously associative. On the other hand, correlation in the frequency domain is multiplication by the complex conjugate, which is not associative.

The associativity of convolution is what allows you to "pre-convolve" the filters, so that you only need to convolve the image with a single filter. For example, let's say you have an image $f$, which you need to convolve with $g$ and then with $h$.  $f * g * h = f * (g * h)$. That means you can convolve $g$ and $h$ first into a single filter, and then convolve $f$ with it. This is useful, if you need to convolve many images with $g$ and $h$. You can pre-compute $k = g * h$, and then reuse $k$ multple times.

So if you are doing template matching, i. e. looking for a single template, correlation is sufficient. But if you need to use multiple filters in succession, and you need to perform this operation on multiple images, it makes sense to convolve the multiple filters into a single filter ahead of time.


Is a real Grayscale image always 8 bit per pixel indexed? I tested several Grayscale images downloaded from the Internet. Some of them are 8bpp indexed while others are 24bpp (true colored). Looks like different Grayscale images have different depths.

For instance, this online application converts a color image to Grayscale. But, when I check its bit depth, I see that bit depth is not changed.

I think, Grayscales with 24bpp are pseudo-grayscales. That is, they are the average of R, G and B components. I think, real Grayscales are always 8 bit per pixel indexed.

Is my notion correct?

If not, please, provide me some links that I can study.

I don't think that it is that easy. Grayscale images do have certain benefits when you think of imaging experiments. When the measurement is squeezed into eight bits, even though the measurement is much more precise, then there is a loss of information due to the quantization. See the DICOM format for example that is widely used in medical imaging. It is very common there that the information is encoded in 12 bit, or even more.


Impulse property - needed mathematical proof Prove:
  $$ \left. \int_{-\infty}^{\infty}\frac{d\delta(t)}{dt}\phi(t)dt=\frac{-d\phi(t)}{dt}\right|_{t=0} $$

$\delta(t)$ is impulse signal and other one is any general signal.

Can anybody do the proof for me?
Or at-least how to proceed further?

For any pair of functions $u(t)$ and $v(t)$ we have that

$$\frac{d}{dt}(u(t)v(t))=v(t)\frac{du(t)}{dt}+u(t)\frac{dv(t)}{dt}$$

If we integrate on both sides of the equation we get

$$\int\frac{d}{dt}(u(t)v(t)) \ dt=\int v(t)\frac{du(t)}{dt} \ dt+\int u(t)\frac{dv(t)}{dt} \ dt$$

In this case, $u(t)=\delta(t)$ and $v(t)=\phi(t)$, so

$$\int\frac{d}{dt}(\delta(t)\phi(t)) \ dt=\int \phi(t)\frac{d\delta(t)}{dt} \ dt+\int \delta(t)\frac{d\phi(t)}{dt} \ dt$$

The left side of the equation can be solved easily as

$$\int\frac{d}{dt}(\delta(t)\phi(t)) \ dt = \frac{d}{dt}\int\delta(t)\phi(t) \ dt=\frac{d}{dt}\phi(0)=0$$

as $\phi(0)$ is a constant. We also know that the delta function has the following property

$$\int \delta(t)\frac{d\phi(t)}{dt} \ dt = \left.\frac{d\phi(t)}{dt}\right|_{t=0}$$

So going back to the first equation and replacing with these results:

$$0=\int \phi(t)\frac{d\delta(t)}{dt} \ dt+\left.\frac{d\phi(t)}{dt}\right|_{t=0}$$

And that's the proof.


Implementation of halfband pass filter I need to interpolate a complex valued bandlimited periodic function using local interpolation.

I can have the signal sampled at any frequency I want over at equispace intervals.

I am aware that for bandlimited periodic functions the best interpolation method is global interpolation (FFT).

However, for this project I need to use local interpolation.
I managed to perform the interpolation using Lagrange polynomials of order 10 but I require an oversampling of 2 to achieve the desired results.
My goal is to reduce oversampling to any extend even if I need a longer interpolator.

I recently read as much as I could about FIR filters but somethings are not clear to me.

Since the interpolation factor is of the order of 2, I am planning to use a halfband pass interpolator (many books say this is the best since one gets a cheap interpolator due to the zeros in the impulse response).

My questions are:

How do I design a halfband pass filter? In Matlab there are nice functions in the dsp package that allow directly to construct such filters, unfortunately I have to do this implementation in C. I searched and I found that the remez() function from Octave to compute the coefficients is available in C and is the one I am planning to use.
Designing a low pass filter with the Parks-Mclellan algorithm is it by default a halfband pass filter (provided I set up correctly the frequency bands) or do I have to enforce the time domain constraints?
Regarding the linear phase of the FIR filter, what I understand is that my signal will be shifted by some samples, how do I compensate to remove the shifting? What comes to my mind is to apply the corresponding shifting on the opposite direction, is this necessary?
If you could provide an working example of an interpolation of a periodic function (e.g. sum of sines) using a FIR I would be very grateful.

The impulse response and frequency response of a length N zero-phase FIR lowpass filter can be designed using the Remez exchange algorithm.

N = 11; % filter length - must be odd

b = [0 0.1 0.2 0.5]*2; % band edges

M = [1 1 0 0 ]; % desired band values

h = remez(N-1,b,M); % Remez multiple exchange design

The impulse response h is returned in linear-phase form, so it must be left-shifted (N-1)/2 = 5 samples to make it zero-phase.
Answer taken from https://www.researchgate.net/post/How_can_I_design_a_null_phase_FIR_filter.

I implemented this solution some time ago and still works very well in my code.


I want to know if there is a practical way to find out if multipath propagation is significant for audio in a given environment? I am employing audio signals generated by a mechanical press for Blind Signal Separation. But, prior to algorithm application, I need to clarify if the signal strengths coming in multipath (except direct path) are significant.

Ideally, I would like to use the four microphone recordings of audio generated by this press for this task. However, if it's difficult or impractical I am interested in trying out other methods.

Cheers!

Never done this but it might work. Record the audio simultaneously using multiple microphones located in the neighborhood of your main recording position. If the cross-correlations between signals from different pairs of microphones show peaks at different time lags than the autocorrelation, then there must be either multiple sources or multiple paths. Be careful not to confuse multiple paths with reflections from the microphone bodies.


Is it impossible to determine the inverse Z-transform without any other information? Suppose I give you this as my transfer function $H(z)$:

$$   H(z) = \frac{1} { 1 - az^{-1}}$$

With no other information given, is it even possible to determine the inverse Z-transform?

The reason I'm asking is because on a particular Z-transform table I have, there are two possible resolutions for this, depending on the desired region of convergence:

As one can see we have two choices,

$a^nu[n]$
$-a^nu[-n-1]$

So I think I need more information about the system to make the right decision, that is, whether the response is supposed to be causal, stable, or both.

It's true that the algebraic expression for the $\mathcal{Z}$-transform is generally not sufficient for computing the corresponding time-domain sequence. The additional information we need is the region of convergence (ROC). The ROC together with the expression for the $\mathcal{Z}$-transform uniquely determines the corresponding time-domain sequence.

For the given example

$$H(z)=\frac{1}{1-az^{-1}}\tag{1}$$

it is clear that the expression $(1)$ can be obtained in two ways. First,

$$\begin{align}-\sum_{n=-\infty}^{-1}a^nz^{-n}&=-\sum_{n=1}^{\infty}a^{-n}z^n\\&=-\frac{\frac{z}{a}}{1-\frac{z}{a}}\\&=\frac{1}{1-az^{-1}},\quad |z|<|a|\end{align}$$

corresponding to the sequence $h[n]=-a^nu[-n-1]$, and, secondly,

$$\sum_{n=0}^{\infty}a^nz^{-n}=\frac{1}{1-az^{-1}},\quad |z|>|a|$$

corresponding to $h[n]=a^nu[n]$. In both cases I've used the formula for the geometric series.

The only difference is the ROC, so it's crucial to know the ROC in addition to the expression for the $\mathcal{Z}$-transform. The ROC $|z|>|a|$ indicates a right-sided sequence, whereas the ROC $|z|<|a|$ corresponds to a left-sided sequence. For higher order systems, you can also have an annulus $|a|<|z|<|b|$ as ROC, corresponding to a two-sided sequence.

If the sequence is interpreted as the impulse response of a system, then the system is stable if the ROC contains the unit circle $|z|=1$. It is causal if the sequence is zero for $n<0$.


I need to insert an underwater acoustic channel model to an existing code that compared between OFDM and GFDM The acoustic channel model that I have will output the channel impulse response. However, my signal is modulated in frequency domain (by doing fft). Hence, I want to obtain a frequency response of the CIR. Question :

Can I perform FFT on the CIR which give me frequency response and then just multiply my signal with it?
If that is correct, is the diagonal of the multiplication will be my signal after passing through the channel?

 H=fft(h,sys.N);    %sys.N = total data symbols h= CIR
 Sig =H.*tf_ofdm;   %tf_ofdm = input signal
 Signal=diag(output);  %output signal

I am doing FFT with the total data symbols because the CIR has random size every time it is being used.

I modeled the the acoustic channel noise in matlab:

    % On the Relationship Between Capacity and Distance in an
    % Underwater Acoustic Communication Channel

    Turbulence=17-(30*log10(f));
    Noise_Turbulence=  power(10,(Turbulence*0.1));

    Shipping=40+(20*(s-0.5))+(26*log10(f))-(60*log10((f)+0.03));
    Noise_Shipping= power(10,(Shipping*0.1));

    Wind=50+(7.5*(power(w,0.5)))+(20*log10(f))-(40*log10((f)+0.4));
    Noise_Wind=power(10,(Wind*0.1));

    Thermal=-15+(20*log10(f));
    Noise_Thermal= power(10,(Thermal*0.1));

    Noise=Noise_Turbulence+Noise_Shipping+Noise_Wind+Noise_Thermal;
Ambient_Noise=10*log10(Noise); %the values will be converted into dB


FFT of a $N$-length real sequence via FFT of a $N/2$-length complex sequence I am doing FFT. The one way I did is to put real part in one array and the imaginary part as zeros and then calculated the FFT. But this take up lots of resources.
Now I want to calculate FFT by Dividing Real Input by 2 in the form of ODD and EVEN samples and put even in real array and odd in Imaginary array. Then calculate the FFT. But after the FFT calculation how can I combine them to get the same results as I got when I used imaginary as zero's.

Any idea??

You can use Decimation In Time (DIT) to calculate the FFT of single $N$ length sequence, using two $N/2$ sequences and combine them later on with a single butterfly.

Knowing that FFT of an even length, real sequence, can be decomposed into even and odd part as:

$$\begin{align}
\mathrm{FFT_N}(x)&=\sum_{m=0}^{N/2-1}x(2m)e^{-2\pi j k(2m)/N} + \sum_{m=0}^{N/2-1}x(2m+1)e^{-2\pi j k(2m+1)/N}\\
&=\sum_{m=0}^{N/2-1}x(2m)e^{\frac{-2 \pi j k m}{N/2}}+e^{-2\pi j k/N} \sum_{m=0}^{N/2-1}x(2m+1)e^{\frac{-2\pi j k m}{N/2}} \\
&= \mathrm{FFT}_{N/2}(x_{even})+e^{-2\pi j k/N}\mathrm{FFT}_{N/2}(x_{odd})
\end{align}$$

and the following holds:

$x_{even}(m)=x(2m)$

$x_{odd}(m)=x(2m+1)$

Now you can form the imaginary signal from both:

$$\hat{x}(m)=x_{even}(m)+jx_{odd}(m)$$
for $m=0\ldots N/2-1$

Calculate the FFT of that complex sequence:

$$\hat{X}=\mathrm{FFT}_{N/2}(\hat{x})$$

Now you can extract corresponding FFT's of both even and odd part:

$$X_{even}(k)=\frac{\hat{X}(k)+\hat{X}(N/2-k)^*}{2}$$
$$X_{odd}(k)=-j\frac{\hat{X}(k)-\hat{X}(N/2-k)^*}{2}$$

And combine them using a butterfly:

$$X(k)=X_{even}(k)+e^{-2\pi j k /N}X_{odd}(k)$$

Keep in mind that those values are for $k=0\ldots N/2-1$. The missing value for $N/2$ can be obtained knowing the symmetry:

$$X(N/2)=\frac{[\hat{X}(0)+\hat{X}^*(0)]+j[\hat{X}(0)-\hat{X}^*(0)]}{2} $$

And that gives you the final result.


Integration of square wave trying to program an integrator. My input is a square wave and my expected output should be a triangle wave. However, whenever I pass it through my low pass filter algorithm (just a 2nd order butterworth low pass filter with a Q of 0.707), I never seem to get a triangle wave. Instead, I get a a smooth square wave (I guess that's similar to a capacitor smoothing it out?). I am not sure how to tackle this problem as I'm a little new to signal processing + algorithms. Any help would be greatly appreciated.

Below is a picture of an example test I did through Xcode. I was sending a F3 note (around 349 hz) through a low pass with a cutoff of 200 hz and the output is shown below. This is probably the sharpest/closest I've gotten to the triangle wave.

EDIT: to clarify, I am trying to convert a square wave to a triangle wave (preferably through a low pass filter if that is totally doable).

If your square wave has a mean of zero (this is important!), then a simple accumulator can do the job. Its operation is described by

$$y[n]=x[n]+y[n-1]$$

where $x[n]$ is the input (square wave) and $y[n]$ is the output (triangular wave).

This is a simple Matlab/Octave script showing how it works:

sq = [1,-1,1,-1,1,-1,1,-1]'*ones(1,5);
sq = sq'(:);                 % a few periods of a square wave
tr = filter(1,[1,-1],sq);    % filtered by accumulator
% plot
subplot(2,1,1),stem(sq)
subplot(2,1,2),stem(tr)


Transfer function, amplitude response and difference equation for a filter I've found a paper with a filter described in terms of transfer function, amplitude response and difference equation:

transfer function of the second-order low-pass filter:

$$
H(z) = \frac{(1-z^{-6})^{2}}{(1-z^{-1})^{2}}
$$

amplitude response (T - sampling period):

$$
|H(\omega T)| = \frac{\sin^{2}(3\omega T)}{\sin^{2}(\omega T/2)}
$$

difference equation of the filter (cut-off 11 Hz, gain 36):

$$
y[nT] = 2y[nT-T] - y[nT-2T] + x[nT] - 2x[nT- 6T] + x[nT- 12T]
$$

What was the process to create these equations and how to run the filter on the signal?

start from the transfer function: $$
H(z) = \frac{(1-z^{-6})^{2}}{(1-z^{-1})^{2}}
$$

expand it: $$
H(z) = \frac{ 1-2z^{-6} + z^{-12} }{ 1 -2z^-1 + z^{-2}  }
$$

Frequency response is:
$$
H(e^{j\omega}) = \frac{ 1-2e^{-j6\omega} + e^{-j12\omega} }{ 1 -2e^{-j\omega} + e^{-j2\omega}  }$$

Note: this may be simplified by proper algebraic grouping and cancellations to what you have provided.

Transform this rational H(z) directly back into its corresponding LCCDE:
$$y[n]-2y[n-1]+y[n-2] = x[n]-2x[n-6]+x[n-12]$$

I can't figure out a meaning for the wT representation, nor if it is necessary...


Slight crackle with overlap-add frequency-domain convolution I'm trying to make an FFT audio filter using overlap-add but there's some slight rhythmic distortion/crackle in the output. I have an irSize (impulse response size) of 512 samples and a blockSize of 512. The overlap is 512 samples. Since fftSize >= blockSize + irSize - 1, fftSize is 1024 samples. I'm using FFTS for the FFT so the output size of forward_real (fftBlockSize) is 1026.

Here are my steps:

// copy input
for(i = 0; i < blockSize; i++)
    fftBuffer[i] = inBuffer[i];

// set rest to zero
for(i = blockSize; i < fftBlockSize; i++)
    fftBuffer[i] = 0;

FFT::forward_real(fftSize, fftBuffer, fftBuffer);

// complex multiplication
for(i = 0; i < irSize + 1; i++)
{
    convBuffer[i * 2] = fftBuffer[i * 2] * irBuffer[i * 2] - fftBuffer[i * 2 + 1] * irBuffer[i * 2 + 1];
    convBuffer[i * 2 + 1] = fftBuffer[i * 2] * irBuffer[i * 2 + 1] + fftBuffer[i * 2 + 1] * irBuffer[i * 2];
}

FFT::inverse_real(fftSize, convBuffer, fftBuffer);

// overlap-add
for(i = 0; i < blockSize; ++i){
    outBuffer[i] = fftBuffer[i] + parts[i];
    parts[i] = fftBuffer[i + blockSize];
}

Here you can see the little jump in the output every 512 samples. This is the source of the distortion.

Any ideas? I've been struggling with this for days. It seems like either I'm doing the complex multiply or the overlap add slightly wrong. But I can't for the life of me figure it out. It seems like I'm doing the theory right.

Thanks.

For FFT fast convolution, the complex multiplication in the frequency domain needs to be done over the full length of the FFT (fftSize), not just the length of the impulse response in the time domain.  Otherwise (for strictly real data and a strictly real impulse response) the IFFT input will likely no longer be conjugate symmetric.


Image resizing and interpolation In digital image processing, I need to understand what does image resizing do? Does it only change the dimensions of the image, for example, from 256×256 to 128×128 with the same pixel spacing (the same dpi)? Or does it change the dpi or pixel spacing in an image? If I need to increase the size of the image and also increase its dpi, how can I do this?

In digital image processing, an image resizing changes only the array size (dimensions) of the image. It does not care about dpi. It just upsamples or downsamples the image array.

DPI (Dot Per Inch) is considered only when the actual physical size of the displayed (or printed) image is of concern.

When you display a scaled image on the same physical monitor, then the apparent dpi of the image content will be changing too. This is especially important on text processing with font rendering.

Some programs try to keep image physical display size the same, while scaling them or when displayed on different monitors, by incorporating a dpi parameter to be used as well.

[edit for answer below]
In those examples from Gonzalez's book, images are considered with their physical sizes when displayed (or printed). Assuming square pixels, DPI of a physical display device is computed by
$$ \text{DPI} = \frac{ \text{ number of horizontal pixels on the monitor}}{ \text{ horizontal length (in inches) of the monitor}} $$

for example for a typical 15.6 inches, FHD (1920 x 1080 ), 16:9 aspect ratio laptop, the DPI will be
$$ \text{DPI} = \frac{ 1920 \times \sqrt{16^2 + 9^2}}{15.6 \times 16}  = 141.212 $$

That's the DPI of the monitor.

Physical size of a displayed image with a array size of $M \times N$ is given by
$$ \text{H} = \frac{ M }{ \text{DPI}} $$ horizontally and

$$ \text{W} = \frac{ N }{ \text{DPI}} $$ vertically in inches...

For amonitor with a given DPI, an image with smaller array size $M,N$ will display smaller on the screen. If you upsample (aka interpolate, or scale) the image array from $M,n$ to, say, $3M, 3N$, then the displayed image size will also be tripled.

And finally; interpolation of an image (or any digital data) means computing the sigal values that's physically in between its existing samples. In that respect you are decreasing the pixel distance when you interpolate. However, interpolation will not give you new data; it will only make exsting data smoother. That's exacty why you cannot get authentic FHD video simply by interpolating from an SD original source...


Plotting the magnitude response of a biquad filter I've put together a multi-band audio EQ using biquad filters.  I'm getting the coefficients using the methods from the RBJ cookbook.

Now I want to plot the curve showing the magnitude response.  I'm using an equation from this source

Here is my function to get the coefficients and get the magnitude response at points of interest.

void GetCoefficients (double samplerate = 44100.0) {

    //from the rbj biquad coefficient cookbook by Robert Bristow-Johnson
    long double SR = (long double)samplerate;
    long double A = powl(10.0L, dBGain/40.0L);
    long double W0 = 2.0L * PI * Center / SR;
    long double alpha = sinl(W0)*sinhl( LN2/2.0L * WidthInOctaves * W0/sinl(W0));

    if (Type == "peaking") {
        b0 = 1.0L + alpha * A;
        b1 = -2.0L * cosl(W0);
        b2 = 1.0L - alpha * A;
        a0 = 1.0L + alpha / A;
        a1 = -2.0L * cosl(W0);
        a2 = 1.0L - alpha / A;
    }

    long double w;
    long double numerator;
    long double denominator;
    long double magnitude;

    for (int i = 0; i < 59; ++ i) {
        w = 2.0L*PI*FreqPoints[i] / SR;
        numerator = b0*b0 + b1*b1 + b2*b2 + 2.0L*(b0*b1 + b1*b2)*cosl(w) + 2.0L*b0*b2*cosl(2.0L*w);
        denominator = 1.0L + a1*a1 + a2*a2 + 2.0L*(a1 + a1*a2)*cosl(w) + 2.0L*a2*cosl(2.0L*w);
        magnitude = sqrtl(numerator / denominator);
        FrequencyResponse[i] = magnitude;
    }
}

My filters sound correct, but my plot seams way wrong.  For example, when I have calculated the coefficients for a peaking filter with a width of $2$ octaves, centered at $4398\textrm{ Hz}$, with $3\textrm{ dB}$ gain; my magnitude response using those coefficients at $349\textrm{ Hz}$ is about $12$.

I must be doing something wrong, but I can't quite figure it out.

Someone over at the KRV forum pointed out that I need to normalize my coefficients so that a0 = 1.

b0 /= a0;
b1 /= a0;
b2 /= a0;
a1 /= a0;
a2 /= a0;
a0 = 1;

Solved my problem.


I am trying to convert a signal into another unit, but the error keeps increasing during each iteration I have a signal, A, and I want to convert it into another unit by using a formula:

$$s_i = \frac{(1+x_i)}{(1-x_i)}\cdot s_{i-1}$$

where $x_i$ is each sample of the signal. I am using a loop where the last value of $s$ is stored and used for the next iteration (as the new value of $P$). It goes on and on until the last sample of the signal.

Signal B is the result of the calculation. However, I noticed that the error keeps increasing during the calculation. As shown in the figure, the red line is how I want the signal to be.

If I understand your notation, you have $y[n] = \frac{1+x[n]}{1-x[n]}y[n-1]$, where $y[n]$ is the output signal.

This can be rewritten as $y[n] = y[n-1] + \frac{2x[n]}{1-x[n]}y[n-1]$. This is an integrator ($y[n] = y[n-1] + w[n]$) where $w[n]$ is some non-linear function of $x$ and $y$.

This integrating component is what gives you the drift in output. Note that the small positive offset in your input signal will result in a drift in your output.

If you could explain your application, and why you came up with the original difference equation, maybe some of us could come up with some insights for a solution.


FSK Demodulation, Bit Timing Recovery I am currently implementing acoustic FSK modulation and demodulation. I am not a signal processing guy so any help about the bit timing recovery would be very appreciated. Currently I  implemented the demodulator using two matched filters for each tone (with a difference of phase of $90^{\circ}$ for non coherent detection). Basically the output of each filter are peaks with different amplitudes. I have two questions:

How can I perform bit timing recovery?
What do you recommend for packet synchronization preamble (chirp,barker code, golden code etc)?

I googled out the Gardner algorithm  but I am not sure if it is applicable and how. Consider I am currently working with 2 frequencies with bitrate of 800 bps using the soundcard.

I am currently implementing acoustic FSK modulation and demodulation. I am not a signal processing guy…

Since you say you have matched filters, and you mention non-coherent detection, I think you're pretty much of a digital communication person already – the step to being a DSP person is pretty small :)
The fully-fledged synchronizer SDR approach
So, the classical way to do this is to take the outputs of your two filters, (magnitude) square them, low pass filter the result to something above $2 f_\text{symbol,max}$, and during/after filtering decimate (if possible), so that we don't waste operations per second.
I made myself a quick real-value bandpass signal FSK demod (I think if I did this again, I would have replaced the two bandpass filters with one combined frequency translation and one low pass filter, giving me complex baseband, or threw a polyphase filterbank at the problem, but whatever) in GNU Radio with the help of the GNU Radio companion (flow graph file here):

Point is that we probably want to symbol timing synchronize between the divide and the threshold block. We could do afterwards – the classical, "microcontroller & digital hardware design" approach would in fact be Muller & Mueller clock sink and that would work there just as well – but let's not waste the slope information.
So, I'd add a polyphase clock sync here; I must admit I could not formulate this better than Tom, a disciple of harris, did in the documentation of the eponymous GNU Radio block:

Timing synchronizer using polyphase filterbanks.
This block performs timing synchronization for PAM signals by minimizing the derivative of the filtered signal, which in turn maximizes the SNR and minimizes ISI.
This approach works by setting up two filterbanks; one filterbank contains the signal's pulse shaping matched filter (such as a root raised cosine filter), where each branch of the filterbank contains a different phase of the filter. The second filterbank contains the derivatives of the filters in the first filterbank. Thinking of this in the time domain, the first filterbank contains filters that have a sinc shape to them. We want to align the output signal to be sampled at exactly the peak of the sinc shape. The derivative of the sinc contains a zero at the maximum point of the $\text{sinc}$ ($\text{sinc}(0) = 1,\, \text{sinc}'(0) = 0$). Furthermore, the region around the zero point is relatively linear. We make use of this fact to generate the error signal.
If the signal out of the derivative filters is $d_i[n]$ for the $i$th filter, and the output of the matched filter is $x_i[n]$, we calculate the error as:
$$e[n] = \frac{\Re\{x_i[n]\} \cdot \Re\{d_i[n]\} + \Im\{x_i[n]\} \cdot \Im\{d_i[n]\}}{2.0}\text{ .}$$
This equation averages the error in the real and imaginary parts. There are two reasons we multiply by the signal itself. First, if the symbol could be positive or negative going, but we want the error term to always tell us to go in the same direction depending on which side of the zero point we are on. The sign of $x_i[n]$ adjusts the error term to do this. Second, the magnitude of $x_i[n]$ scales the error term depending on the symbol's amplitude, so larger signals give us a stronger error term because we have more confidence in that symbol's value. Using the magnitude of $x_i[n]$ instead of just the sign is especially good for signals with low SNR.
The error signal, $e[n]$, gives us a value proportional to how far away from the zero point we are in the derivative signal. We want to drive this value to zero, so we set up a second order loop. We have two variables for this loop; $d_k$ is the filter number in the filterbank we are on and $d_\text{rate}$ is the rate which we travel through the filters in the steady state. That is, due to the natural clock differences between the transmitter and receiver, $d_\text{rate}$ represents that difference and would traverse the filter phase paths to keep the receiver locked. Thinking of this as a second-order PLL, the $d_\text{rate}$ is the frequency and $d_k$ is the phase. So we update $d_\text{rate}$ and $d_k$ using the standard loop equations based on two error signals, $d_\alpha$ and $d_\beta$. We have these two values set based on each other for a critically damped system, so in the block constructor, we just ask for "gain," which is $d_\alpha$ while $d_\beta$ is equal to $\frac{\text{gain}^2}{4}$.

Realizing this is audio frequency
So, parameterizing this block properly (and probably reducing bpf_decim andlpf_decim to give the synchronizer more headroom to shift phases) you could build a very stable symbol timing recovery, that would be totally overkill for your application
Since you're doing sampling with a sound card on both ends, and since I presume you don't have to deal with Doppler, the symbol rate has a fixed relation to the carrier frequencies (e.g. the transmitting side has a higher carrier with a period of 8 sound card samples, and a symbol always takes let's say 128 samples, so there's a fixed ratio), you can take a short cut:
If you know the ratio between carrier frequency and symbol rate, just estimate one and you get the other for free! In this case, FSK, that's particularly easy: just wait for a period when the e.g. upper band is active, and use a quadrature demodulator (i.e. practically $\frac{d\arctan}{dt}[n]$) to give you an estimate of the carrier frequency after your decimating band pass filter, and use the resulting value as frequency correction factor – you immediately solve the carrier frequency recovery and symbol rate recovery problem at once. I don't even think you'd need timing recovery after that – after a thresholding, if you oversample enough (which, at these rates, isn't a problem), you can just make a majority decision in a window, or shift that window as long as necessary to get a maximum in the window.
Regarding Packet Sync
Well, since you've got a working demodulator that will probably run too much out of sync when there's no signal present, just use a fixed known data sequence to find your packet.


How to interpret frequency response of moving average system in discrete time signal processing? Excerpt from Discrete Time Signal Processing by Alan Oppenheim:

The text says this can be considered as an approximation of a Low Pass Filter. A low pass filter in continuous signal system is something that would block frequencies above a cutoff. But here in the discrete system, how is this to be understood? In a discrete time system, "frequency" is just angle/sample. Does the picture above mean "the system allows samples taken at the rate (0 to $2\pi/5$) / sample and blocks the rest" ?
EDIT:
Is it like this: say there is signal with a frequency range 20Hz to 50Hz. The sampling frequency in 100Hz. The lower frequency components get smaller $\omega$ and for higher frequency (50Hz) it gets close to $\pi$. If I have a low pass filter, it would block out those samples with large $\omega$. And this functionality now resembles a continuous time low pass filter. Is this a correct interpretation?

In discrete time, a frequency of zero is just a constant sequence, in complete analogy with continuous time. The main difference between discrete time and continuous time is that in discrete time you have a maximum frequency, whereas in continuous time you (theoretically) don't have one. That maximum frequency in discrete time is achieved for a sequence with alternating values, i.e., the sequence $(-1)^n$ has a maximum frequency of $\omega=\pi$. There can be no larger frequency.
If the sequence was obtained from sampling a continuous-time signal with a sampling frequency $f_s$, then $\omega=\pi$ translates to an actual frequency in Hertz via $f=\omega/(2\pi f_s)$.
Returning to the example of the discrete-time lowpass filter, from the frequency response you see that higher frequencies (up to the maximum frequency $\omega=\pi$) are attenuated. That means that a discrete-time sinusoid's amplitude is changed according to the magnitude of the filter's frequency response at the sinusoid's frequency.
So if the input signal equals
$$x[n]=A\sin(n\omega_0)\tag{1}$$
the amplitude of the corresponding output signal is given by $A|H(e^{j\omega_0})|$. If $\omega_0$ happens to be a multiple of $2\pi/5$ then the input is completely suppressed because these frequencies are the zeros of the filter's frequency response.


what are frequencies in a digital image? I don't completely understand the term " frequency" with respect to digital image processing.
for 1D electrical signals , its pretty easy to understand , visually too.
higher the frequency - more packed( or closer) the signal looks in a given time interval than its lower frequency counterpart.

But when it comes to images , with parameters like intensity, 8-bit gray scale, pixel value,etc. - what does frequency mean ?

In images, intensity means amplitude and you can see it as in 1D case.
More packed (or change of intensity) in 2D area is higher frequency.
Note that, this way is just visual approximation by human's eyes.


Does it make sense to express a given FT from 0 to 1 as a convolution? I'm learning about Fourier transform, and was asked whether the following FT can be expressed as a convolution: $$X[k] = \sum\limits_{n=0}^{N=1} x[n]e^{-i\frac{2\pi}{N}kn}$$
There are two things I don't understand:

Why would anyone want to express the FT itself as a convolution? from what I understand, FT is used to express convolution, as stated by the Convolution Theorem.
Assuming it does make sense to express an FT as a convolution, how can we do that with the above FT?

As for the why, it's always good to come up with alternatives, even if one doesn't immediately see their benefit. In this case it's about efficiency, especially if the FFT length is a prime. The resulting algorithm can also be used to solve the more general problem of computing (samples of) the $\mathcal{Z}$-transform on circles or spirals in the complex plane. You can also zoom into a certain frequency region without the need to compute the DFT over the whole range.
The resulting algorithm is called Bluestein's algorithm or Chirp Z-transform. The basic trick to represent the DFT (or $\mathcal{Z}$-transform) by a convolution is to write
$$nk=\frac{n^2+k^2-(n-k)^2}{2}\tag{1}$$
and use this expression in the original formula. This results in the representation of the original DFT sum by $(i)$ a multiplication of the original sequence with a chirp sequence, $(ii)$ a convolution, and $(iii)$ another multiplication with a chirp sequence. The details are explained in the link cited above.


Image processing (number recognition) woes So I have the following image:

Using GIMP, I was able to get to here:

Select all white color (with a touch of fuzz)
Invert selection
Cut selection
Invert all color (goes black)

But now I'm stuck...

I want to use something like ocrad to extract the numbers.

ocrad goes crazy when I input the black&white image above.

I'd like to see (the grammar is the numbers 1 to 45 inclusive):

var a = [30,31,32,35,37,40,44];
var b = [6,7,11,15,18,21,22];
var c = [5,11,15,18,23,37,28];

It's worth pointing out that the image

is always the same dimensions
is always the same orientation
therefore, I think I can split the image into three rows relatively robustly

Like so:

(in the above case, I also did an extra step - remove black contiguous regions that have large area. Though I don't know how to do this programatically - I did this by eye in GIMP using the "magic wand" tool).

The above image does OK in ocrad, but it gets into trouble with the ring surrounding the "44".

I'm trying to avoid pattern matching.

Does anyone have any tips? I'm thinking of doing all this in Python/OpenCV/similar.

You may want to check PyTesser. Also, this article is worth reading.


Clarification for "Digital Signal Processing" by Steven W. Smith related to complex representation of sinusoids I have a question related to complex number representation of sinusoids from a book by Steven W. Smith "Digital Signal Processing". The book is made freely available by the author on dspguide.com. The part I have question about is available here, pp.559-563.

Excerpt from the book (p. 560):
A sinusoid is represented by a complex number in rectangular form as follows:
$$
A \cos(\omega t) + B \sin(\omega t) \rightleftharpoons a + jb
$$
where $a \rightleftharpoons A$ and $b \rightleftharpoons -B$.
A sinusoid is represented by a complex number in polar form as follows:
$$
M \cos(\omega t + \phi) \rightleftharpoons M e^{j\theta}
$$
where $\theta \rightleftharpoons -\phi$.
Why change the sign of the imaginary part and phase angle? This is to make the substitution appear in the same form as the complex Fourier transform described in the next chapter. The substitution techniques of this chapter gain nothing from this sign change, but it is almost always done to keep things consistent with the more advanced methods.

What I am confused about is that these two representations are not equivalent. Here is an example following the rules described above for rectangular and polar forms:
$$
\cos(\omega t) + \sin(\omega t) \rightleftharpoons 1 - j = \sqrt{2} e^{-j\frac{\pi}{4}}
$$
$$
\cos(\omega t) + \sin(\omega t) = \sqrt{2} \cos(\omega t - \frac{\pi}{4}) \rightleftharpoons \sqrt{2} e^{j\frac{\pi}{4}} = 1 + j
$$
The rules for rectangular and polar form produce phase angles of opposite sign.
Author gives few more examples on p. 561 (first paragraph in "Complex Represenation of Systems") that show clear discrepancy between rectangular and polar forms considering phase angle.
I am aware that we can take either rectangular or polar form to "be correct", the sign is a convention anyway. But how can both of them be used at the same time?

Phasors are useful for the analysis of (real-valued) linear time-invariant (LTI) systems. A phasor is a complex number
$$C=|C|e^{j\phi}=A+jB\tag{1}$$
which represents a sinusoidal signal:
$$x(t)=\textrm{Re}\big\{Ce^{j\omega t}\big\}=A\cos(\omega t)-B\sin(\omega t)=|C|\cos(\omega t+\phi)\tag{2}$$
That's already the end of the story. Eq. $(2)$ is the generally accepted definition of a phasor. The definition in the book is either a mistake or an idiosyncrasy making things unnecessarily confusing.


Digital Signal Processing DSP Overview study documents Well, I couldn't find better expression to say what I want :(

I'm looking for a DSP study source, that doesn't dive too deep, and covers "all" or most of DSP topics, I just want to have an overview of it, then, as a next step, it will be easier to go deeper.

A PDF file, an online file, video course, tutorial, ... I french we say "Résumé de cours" or "course summary" but it didn't give me what I expected!

Plz, for those having no answer, plz, tell me how to ask about it and google it

Thanks in Advance

There are some books (some free) that cover most DSP topics with a claim of easy approach, probably they cover the most important fundamentals:

1- Understanding Digital Signal Processing_R.Lyons

2- Scientist and Engineer's guide to DSP_S.W.Smith (free from web)

3- Signal Processing First_McClellan


Hide text in audio and got text by phone microphone when the audio is plaid anywhere Thanks for your time to read my question first. I’m new about audio, but I have to make this done:

Add a few bytes of data to a sound file and generate a new sound file.
If I convert that file to any sound format, these bytes data should still there.
When I play that file anywhere, I pick up the sound with mobile phone microphone and got these bytes back.

I have searched some informations:

It's like "Audio Watermark"
A study : HIDING TEXT IN AUDIO USING MULTIPLE LSB STEGANOGRAPHY AND PROVIDE SECURITY USING CRYPTOGRAPHY
This video is what I want, but I'm still not clear how to do that. Could you please give me some suggestion about that?

Thank you in advance for any help you might be able to provide.

I'm not sure about how much experience you have in writting software, but if what you want is to hide some data information within an analog audio signal, it seems for me like the way to go would be to use spread spectrum modulation and encode the data digitally and add it to your audio signal as some background noise.
Spread spectrum is in itself already a form of encryption, so the bytes you add there will be from the get go, very difficult to recover for a hacker unless they: 1-know that are there, and 2-have the key for generating the pseudorandom sequence to de-spread the signal.
But if you want even more security you can choose to packet your data by SSH tcp packets, and have in your phone a program to extract those from the digital signal.

As for point number 2. My suggestion immplies that you add something to the analog signal, or if you will to the PCM raw waveform, which is the same for all practical purposes. So, taking into account that the job of an encoder/compressor like mp3 is to compres the signal, but to eventually recover the original waveform with as less RMS error as possible, we could say that the original noise waveform that you add will still be there, slightly modified, but still there, and hence when you recover the bits by demodulation, and apply suitable Error correction tecniques, I think your data should still be there intact as it was put into the sound.

Please take into account that i just suggested a possible technology, normally used with narrowband RF signals, but maybe you could do some research and check if it will work in your scenario.

Also Another technique that could be used, if you want to send Text (Like they do in the video), is to generate based on your text a signal that will look like the ASCII characters in the time-frequency plane. Then to recover you just perform time-frequency analisys and use any image recognition algorithm.
This just came to mind, so it probably sounds very wierd, but may be another way to add another signal to an audio file.

Hope this helps.


Expression to represent alternating pulse train? I have a question regarding an alternating pulse train signal, $p(t)$, and a system which looks like this:

Immediately, I'm supposed to be able to extrapolate that
$$
P_t = P_1(t)-P_1(t-\Delta)
$$

How do you obtain this?
How would you get the Fourier transform of this pulse train from there?

Since this is a periodic signal, we have to first represent the signal $x(t)$ in terms of its Fourier Series and then compute its Fourier transform

$$\text{Representing } x(t) \text{ in terms of its Fourier series }$$

Let us assume $\Delta=T$

$$T_0=2T$$

$$x(t)=\sum_{n=-\infty}^\infty  x_ne^{j2\pi \frac{n}{T_0}t}=\sum_{n=-\infty}^\infty  x_ne^{j\pi nt}$$

$$x_n=\frac{1}{T_0}\int_{\alpha}^{\alpha +T_0} x(t)e^{-j2\pi \frac{n}{T_0}t}dt=\frac{1}{2T}\int_{0T^-}^{2T^-} x(t)e^{-j2\pi \frac{n}{2T}t}dt$$

$$\frac{1}{2T}\int_{0T^-}^{T^-} \delta (t)e^{-j\pi \frac{n}{T}t}dt+\frac{1}{2T}\int_{1T^-}^{2T^-} -\delta (t-T)e^{-j\pi \frac{n}{T}t}dt$$

$$=\frac{1}{2T}(1)-\frac{1}{2T}e^{-j\pi \frac{n}{T}T}=\frac{1}{2T}(1-e^{-j\pi n})=\frac{1}{2T}(1-(-1)^n)$$
$$x(t)=\sum_{n=-\infty}^\infty \frac{1}{2T}(1-(-1)^n)e^{j\pi\frac{n}{T}t} \tag 1$$

$$\text{Determing Fourier Transform of } x(t) \text{ from its Fourier series representation}$$

Let $n=2k+1$, then (1) becomes
$$=> x(t)=\sum_{k=-\infty}^\infty \frac{1}{2T}(2)e^{j\frac{(2k+1)}{T}\pi t}$$

$$X(jw)=\int_{-\infty}^\infty x(t)e^{-jwt}dt$$
$$=\int_{-\infty}^\infty \sum_{k=-\infty}^\infty \frac{1}{2T}(2)e^{j\frac{(2k+1)}{T}\pi t}e^{-jwt}dt$$

$$\frac{1}{T}\sum_{k=-\infty}^\infty \int_{-\infty}^\infty [1.e^{j\frac{(2k+1)}{T}\pi t}]e^{-jwt}dt$$

$$\text{ Since } F[1]=2\pi\delta (w) \text{ and } F[e^{j\Omega_0t}]=2\pi \delta(\omega-\Omega_0) \text{ Therefore}$$

$$X(jw)=\frac{1}{T}\sum_{k=-\infty}^\infty 2\pi \delta(w-\frac{(2k+1)}{T}\pi)$$

$$=\frac{2\pi}{T}\sum_{k=-\infty}^\infty  \delta(w-\frac{(2k+1)}{T}\pi) \tag 2$$

$$X(f)=\frac{2\pi}{T}\sum_{k=-\infty}^\infty  \delta(2\pi f-\frac{2(2k+1)}{2T}\pi)$$

$$\text{Since } \delta(ax)=\frac{1}{|a|}\delta(x)$$

$$X(f)=\frac{2\pi}{|2\pi|T}\sum_{k=-\infty}^\infty  \delta( f-\frac{(2k+1)}{2T})$$

$$=\frac{1}{T}\sum_{k=-\infty}^\infty  \delta( f-\frac{(2k+1)}{2T})$$
Let $2k+1=m$
$$X(f)=\frac{1}{2T}\sum_{m=-\infty}^\infty (1-(-1)^m)\delta(f-\frac{m}{2T}) \tag 3$$


Extending 1D window functions to 3D (or higher) (cross post of SO: https://stackoverflow.com/questions/27345861/extending-1d-function-across-3-dimensions-for-data-windowing )

For the sake of image (volume) registration, I'd like to apply a windowing function to input data, such that the non-periodic image boundaries do not cause streaking in the FFT. I'm using the example from here for 2D data:

http://mail.scipy.org/pipermail/numpy-discussion/2008-July/036112.html

h = scipy.signal.hamming(n)
ham2d = sqrt(outer(h,h))

Is this extensible to 3D or even N-D?

this question is old, but  I needed same function.

I wrote a function that receive one-dimensional window as parameter and transformed into a cube window

function [win] = window3D(w)

    L = size(w,1);
    m1 = w(:)*w(:).';
    win1 = repmat(m1, [1 1 L]);
    m2 = w(:)*ones(1,L);
    win2 = repmat(m2, [1 1 L]);
    win2 = permute(win2, [3, 2, 1]);
    win = win1.*win2;
end


Power Amplifier Characterization I am trying to characterize a power amplifier fitted to transmitter board.
I applied a sinewave at the required carrier frequency and I collected the input samples and the output samples by using a digital scope. I have noticed a strange behaviour which I can not explain.

when the input sample is zero the output is not zero but a significant value.

Should I expect when the input sample is zero, the output sample from the power amplifier should be zero ?

May the input and the output samples are not in sync.
I tried to sync them, but still the zero input samples do not correspond to very small signal or zero sample.

When I plotted the output/input  I received this plot.

Is this an expected plot ?

How I can get a smooth plot for the ratio of the output over the input ?

I appreciate any comments.

Yes, this is (kind of) expected:

No semiconductor device in itself is truly linear over its full damage-free range; nor is it memoryless.

You're very likely seeing spikes whenever your input approaches zero, so that your division becomes unstable (mathematically), considering you'd still be amplifying noise and an offset current.

Also, you're likely seeing sub-linear amplification when your instantaneous signal becomes large. For a proper discussion I'd just refer you to a research on "intermodulation products", which you'll find in every real (non-perfect) amplifier.

To give you a short glimpse: Let's assume your amplifier is basically a transistor with support circuitry.

A transistor's output current (assume that's emitter current) over input voltage (assume that's base-emitter voltage) curve is (simplifiedly, not considering temporal effects etc):

$$I_\text{E} = I_\text{ES} \left(e^\frac{V_\text{BE}}{V_\text{T}} - 1\right) \tag1\label 1$$

Where $V_\text{T}$ and $I_\text{ES}$ are temperature- and semiconductor dictated constants.

You'll notice this is nothing close to a linear function. Now, the job of an amplifier designer is to operate the components (mainly, transistors) in the amplifier in a way that looks as linear as possible (for power amplifiers, that is a competing requirement to efficiency and large power capability, so these aren't typically all to linear or low-noise).

Imagine you do a Taylor series expansion of $\eqref 1$. You'll notice that the exponential function has the series representation

$$e^x = \sum\limits_{n = 0}^{\infty} \frac{x^n}{n!}\text.\tag2\label2$$

Transistor amplifiers aren't really perfectly represented by that exponential function, as you can do things like influence "weights" on individual series elements within bounds by careful semiconductor design:

$$f(x)=\sum\limits_{n = 0}^{\infty} \tilde a_n \frac{x^n}{n!}=\sum\limits_{n = 0}^{\infty} a_n {x^n}\label 3$$

What you don't want is the $x^0$ term, because that's just a constant bias current ($x^0 \equiv 1$), so you'd typically try to have a very small $a_0$.

What you do want is the $x^1=x$ term, because that's your linear amplification – and typically, you'd want that to be large, so if you can, you see towards that $a_1$ becomes large.

$a_2$ and higher lead to intermodulation, and you want them to be small – but compensating them technically leads to limited $a_1$, and that means your amplifier doesn't amplify much.

Now, your amplifier manufacturer won't offer you a specific input-output curve, often, but other things more important for people that work their amplifier with signals of different frequencies:

Second and Third order intercept points. I can only recommend reading these articles closely, and possibly get a book on semiconductor practice (Tietze-Schenk would be the go-to bible where I'm from), as you're working with amplifiers and need to know what you're doing


Matlab produces two unknown spikes in custom FFT I'm making a relatively small Cooley–Tukey FFT in Matlab and I'm noticing unusual spikes in the result compared with Matlab's own FFT.

The figure below shows the signal flow of my program. It's a standard Cooley-Tukey scheme.

My results when computing the FFT of a 16 Hz sinusoid are shown below

And Matlab's own FFT is shown below

Clearly I'm missing something. What could be the cause of the two large spikes in the middle? I conjecture that it has something to do with how the even and odd parts are combined, for instance if there is a discontinuity there. I'm really not sure though.

Any help is greatly appreciated!

The code I'm using is as follows

clear all

% Generate input data sequence and plot
N=128;
f1=16;
num_cycles=2;
fs=f1*N/num_cycles;
x_time=0:1/fs:num_cycles/f1-1/fs;
x=sin(x_time*2*pi*f1);
plot(x_time,x);

% split inputs into even and odd samples and compute fft of each division
X_o=x(1:2:N);
X_e=x(2:2:N);

fft_x_o=fft(X_o);
fft_x_e=fft(X_e);

% Generate base twiddle factor
W32=exp(-1i*2*pi/32);

% Combine fft even and odd with twiddle factors to produce final output
for k=0:N-1
    if k<N/2
        X(k+1)=fft_x_e(k+1)+(W32^k)*fft_x_o(k+1);
    else
        X(k+1)=fft_x_e(k+1-N/2)+(W32^k)*fft_x_o(k+1-N/2);
    end
end

% plot butterfly fft and matlab fft
FFT_xaxis=0:fs/N:fs-fs/N;
figure
plot(FFT_xaxis,abs(X))
title('Butterfly FFT')
xlabel('Frequency')
ylabel('Magnitude')
matlab_fft=fft(x);
figure
plot(FFT_xaxis,abs(matlab_fft))
title('Matlab FFT')
xlabel('Frequency')
ylabel('Magnitude')

The following is the corrected code. It seems you have the problem in the twiddle factor and the selection of even and odd samples of x[n]...

N=128;
f1=16;
num_cycles=32;
fs=f1*N/num_cycles;
x_time=0:1/fs:num_cycles/f1-1/fs;
x=sin(x_time*2*pi*f1);
plot(x_time,x);

% split inputs into even and odd samples and compute fft of each division
X_o=x(2:2:N);            % odd samples begin at x(2) --> x[1] in sequence
X_e=x(1:2:N);            % even samples begin at x(1)  --> x[0] in sequence

fft_x_o = fft(X_o, N/2);
fft_x_e = fft(X_e, N/2);

% Generate base twiddle factor
WN=exp(-1i*2*pi/N);

% Combine fft even and odd with twiddle factors to produce final output
for k=0:N-1
    if k<N/2
        X(k+1) = fft_x_e(k+1)+(WN^k)*fft_x_o(k+1);
    else
        X(k+1) = fft_x_e(k+1-N/2)+(W3N^k)*fft_x_o(k+1-N/2);
    end
end

% plot butterfly fft and matlab fft
FFT_xaxis=0:fs/N:fs-fs/N;
figure
plot(FFT_xaxis,abs(X))
title('Butterfly FFT')
xlabel('Frequency')
ylabel('Magnitude')
matlab_fft=fft(x);
figure
plot(FFT_xaxis,abs(matlab_fft))
title('Matlab FFT')
xlabel('Frequency')
ylabel('Magnitude')


finding frequency range and dominant frequency of a stored signal I have sEMG signals, captured from commercially available device. the signals are stored in .tdms file format. i want to varify whether the captured signals are really sEMG (or simply EMG) signals or not?
typical frequency range of EMG signals is 5 to 500 Hz and dominant frequency is 50 -150 Hz. can anyone guide me how may i find this from a recorded signal, using either MATLAB or LabVIEW.

To detect the fundamental frequency of your signals, you can use self-correlation. The maximum of the self-correlation will correspond to your fundamental frequency (the second maximum, as self-correlation is always maximum at 0). You can then interpret the lag value of this secondary maximum to find the periodicity of your signal, and the amplitude of the maximum to evaluate the strength of periodicity.
note : The value of the self-correlation at lag=0 is the energy of your signal. This can help as a scale for the periodicity strength.

Otherwise, to visually analyse the spectral content of your signal, you can just apply a Fourier Transform to it and plot the magnitude of the Fourier Transform.


Can I simplify $x=\frac{\ln(|\mathcal F((\mathcal F^{-1}(\sqrt{(f^{-5/3})}))^{2})|)}{\ln(f)}$ Given a variable in time $u_{k}$, and an other variable f which represents the frequencies in a range [a,b]:
$$|\mathcal F \left \lbrace u_{k}\right \rbrace |^{2}=f^{-5/3}\tag{1}$$
where $\mathcal F$ stands for fast Fourier transform.
Now, I want to find $x$ knowing that:
$$|\mathcal F \left \lbrace u_{k}^{2}\right \rbrace|= f^{x}\tag{2}$$
I can isolate $x$ as:
$$x=\frac{\ln(|\mathcal F \left \lbrace u_{k}^{2}\right \rbrace |)}{\ln(f)}\tag{3}$$
I can find $u$ using (1):
$$u=\mathcal F^{-1}(\sqrt{(f^{-5/3})}) \tag{4}$$
On the consequence $x$ is:
$$x=\frac{\ln(|\mathcal F \left \lbrace ( \mathcal F^{-1} \left \lbrace \sqrt{(f^{-5/3})} \right \rbrace )^{2} \right \rbrace|)}{\ln(f)} \tag{5}$$
Is there any way in which I can simplify $x$?? I am stuck with this problem for a lot of weeks...

Here is a potential starting point: it is considered  that when a function behaves as a power law  ($x(t) \sim t^\alpha, t> 0$), then its magnitude spectrum as well, with $|X(f)| \sim f^{-(1+\alpha)}$. One can see for instance:

Computing Fourier transform of power law
An Interesting Fourier Transform - 1/f Noise
Sketches of a "proof" and references (since there are  convergence issues) are provided in an answer to I cannot find reference (paper) of this relation $u(t)t^{α} ↔^{FT} f^{-(α+1)}$.

If I plug that in your case, since $|X(f)|^2 \sim f^{-5/3}$, so we seek $-5/6=-(1+\alpha)$. I obtain $\alpha = -\frac{1}{6}$. So you would have $x^2(t) \sim t^{2\alpha}$, with a magnitude spectrum $|\chi(f)|$ in $f^{-(1-\frac{1}{3})}$, or $f^{-\frac{2}{3}}$.
Now, is the discrete signal sufficiently well-sampled and clean to verify it in the discrete domain is another question.


Modulation with Maximum Likelihood decoder and capacity My question comes from the high spectrum efficiency claim of OFDM. This seems not correct IMHO, at least theoretically.

For me, we lose energy in using cyclic prefix, hence OFDM should be a waste of theoretically achievable rate. That said, OFDM has advantage in implementation and it is easier to get a high throughput.

Long story short, my question is: if complexity is not my concern, can (any) modulation (without channel coding) combined with a Maximum Likelihood (or Maximum a Posteriori) decoder achieve asymptotically the Shannon capacity ?

You are correct that OFDM has several drawbacks, such as the required guard band, number of pilots, and its high PAPR. Many alternatives have been developed, with different sets of pros and cons. Implementation simplicity remains OFDM's most important advantage.

One thing to keep in mind is that OFDM is designed for the wideband (frequency selective) fading channel. The capacity of such a channel is not the usual $C=B\log(1+SNR)$, which only applies in the AWGN channel. Since the fading channel is time-varying, other measures of capacity (such as outage probability) become more important.

To answer you specific question: to my knowledge, it is impossible to reach capacity without coding, except when $\text{SNR} \rightarrow \infty$.


Error implementing Harmonic Product Spectrum algorithm I'm trying to find the fundamental frequencies present in a piano recording using MATLAB. These are the steps I've followed;

Find the envelop of the signal
Find the note onsets
perform FFT between each onset
Harmonic product spectrum.

It's when I try to implement the HPS algorithm that I face a "dimensions don't agree" error. This is the whole code that implements the HPS algorithm .

**% In harmonic prodcut spectrum, you downsample the fft data several times and multiply all those with the original fft data to get the maximum peak.
    %HPS
    seg_fft = seg_fft(1 : size(seg_fft,1)/2 );
    seg_fft = abs(seg_fft);

%HPS: downsampling
for i = 1:length(seg_fft)
    seg_fft2(i,1) = 1;
    seg_fft3(i,1) = 1;
    seg_fft4(i,1) = 1;
%   seg_fft5(i,1) = 1;
end
for i = 1:floor((length(seg_fft)-1)/2)
    seg_fft2(i,1) = (seg_fft(2*i,1) + seg_fft((2*i)+1,1))/2;
end
for i = 1:floor((length(seg_fft)-2)/3)
    seg_fft3(i,1) = (seg_fft(3*i,1) + seg_fft((3*i)+1,1) + seg_fft((3*i)+2,1))/3;
end

for i = 1:floor((length(seg_fft)-3)/4)
    seg_fft4(i,1) = (seg_fft(4*i,1) + seg_fft((4*i)+1,1) + seg_fft((4*i)+2,1) + seg_fft((4*i)+3,1))/4;
end

%HPS, PartII: calculate product
p1 = (seg_fft3)  .* (seg_fft4);
p2 = p1.* (seg_fft2);
p3 = p2.* (seg_fft);

HPS, PartIII: find max
[f_y1,I] = max(p3)

 for c = 1 : length(p3)
     if(p3(c,1) == f_y1)
         index = c;
     end
 end

 % Convert that to a frequency
 f_y(h) = (index / NFFT) * FS;

 h = 1;
 h=h+1;
f_y = abs(f_y)';

 end**

Before implementing p3 = p2.* (seg_fft); the sizes of seg_fft, seg_fft2, seg_fft3, seg_fft4 all have the same dimensions of 16384 1. Then when I DO implement p3 = p2.* (seg_fft); the size of seg_fft changes to 8192 1 while the sizes of rest remain at 16384 1 thus causing an error in multiplication as the dimensions aren't the same.

I'm really confused as to why this keeps happening and nothing I try seems to work out. Would really REALLY appreciate some help here... If someone could fix this code it'd be a GREAT help... Thanx in advance.. I'm real desperate here......

Pad the high frequency portions of all the down-sampled vectors with either zeros, or some noise floor estimate, depending on your desired frequency search range, to make them all the same length or dimension.


Deriving the frequency spectrum graph from the graph of a spectrogram and vice versa? I have a signal's spectrogram and frequency spectrum (which, if i'm not mistaken, is its the signal's fourier transform) as given below:

I understand that the spectrogram plots the frequency (y-axis) vs time (x-axis) and that the colors represent the amplitude.  I also understand the frequency spectrum as the plot of amplitude vs. frequency.  But I'm not quite sure how to produce one graph from the other.   Since there is no time axis on the frequency spectrum, surmise that it is impossible to create the spectrogram from a frequency spectrum.  It seems plausible to derive the frequency spectrum from the spectrogram, but I'm not quite sure how to remove the time parameter.
Any help with this would be greatly appreciated.

The data that when colored produces one vertical slice of your spectrogram image is a spectrum.  If you have the color mapping function, and it is invertible, you could take the array of colors from one vertical column in a spectrogram and attempt to produce the original spectrum for that point in time.  But the inversion process would likely be lossy, since there are usually less quantized colors than spectrum amplitudes.

Since a spectrum is only one vertical slice of a spectrogram, a spectrogram produced from just one spectrum would only be something like 1 pixel wide, thus not very useful.  But a bunch of spectrums at evenly spaced time increments could be used to produce a wider spectrogram.

Also note that due to time-frequency-resolution trade-offs, the width of time data used for each FFT/DFT spectrum array can be different between spectrum and spectrogram analysis presentations.


How can I detect pipeline cracks using OpenCV and Python? I have developed a robot that captures images of the pipeline interior as it moves. The requirement was to be able to detect cracks inside. So far i tried several OpenCV codes that find the crack contours but i was not successful.

Code I'm working on:

import cv2
import numpy as np
image = cv2.imread('pipe_photo1.jpg')
blurred = cv2.pyrMeanShiftFiltering(image,41,91)
gray = cv2.cvtColor(blurred,cv2.COLOR_BGR2GRAY)
ret, threshold = cv2.threshold(gray,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
_, contours, _ = cv2.findContours(threshold, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)
print (len(contours))

cv2.drawContours(image,contours, -1,(0,0,255),6)
cv2.namedWindow("Display",cv2.WINDOW_NORMAL)
cv2.imshow("Display",image)
cv2.waitKey()

Below is the image i obtained from the camera. I want to detect only the crack shown at the bottom of the pipe and be able to draw it using red lines. Your help will really save me in achieving my objectives before its due.

Try this:

Don't convert the image to grayscale, instead convert to HSV
Then, by changing the range of H,S and V, try to mask the crack (using a trackbar can be very helpful)

For example:

lower_value_blue = np.array([80,106,5])    # creating a range of HSV colours
higher_value_blue = np.array([150,255,150])
mask1 = cv2.inRange(blur1,lower_value_blue,higher_value_blue)

This one mask blue colour, without using trackbar

Hope you will get the result


What does "two signals aliasing to each other" mean? So I just read something that two signals are aliasing to each other. I know what a single signal is being aliased - that is when sampling rate is not higher than or equal to the Nyquist rate. But I am just confused of what two signals aliasing to each other mean.

Let's say that each signal only carries single frequency. Let $f_s$ be sampling frequency and $f_1$ be the first signal and $f_2$. $f_s$ is assumed to be lower than Nyquist rate for both signals. In such case, in what condition of $f_1$ and $f_2$ would aliasing occur, given $f_s$?

Aliasing would occur if $f_2 = f_1 + kf_s$, where $k$ is an integer. Consider the following example:

$$
x_1(t) = \cos(j2\pi f_1 t + \phi_1) \\
x_2(t) = \cos(j2\pi f_2 t + \phi_2)
$$

An ideal sampler samples the above at a sample rate of $f_s$ to yield discrete-time signals $x_1[n]$ and $x_2[n]$:

$$
x_1[n] = \cos\left(2\pi f_1 \frac{n}{f_s} + \phi_1\right) \\
x_2[n] = \cos\left(2\pi f_2 \frac{n}{f_s} + \phi_2\right)
$$

We got the above by noting that the sampler takes samples at times $\left[0, \frac{1}{f_s}, \frac{2}{f_s}, \ldots\right]$. To see the aliasing effect, substitute back in the relationship between $f_1$ and $f_2$ that we started with:.

$$
x_2[n] = \cos\left(2\pi (f_1 + kf_s) \frac{n}{f_s}+ \phi_2\right) \\
x_2[n] = \cos\left(2\pi f_1 \frac{n}{f_s} + 2\pi kf_s \frac{n}{f_s}+ \phi_2\right) \\
x_2[n] = \cos\left(2\pi f_1 \frac{n}{f_s} + 2\pi kn + \phi_2\right)
$$

Now we recall a property of sinusoids:

$$
\cos(x + 2\pi m) = \cos(x)\ \forall m \in \mathbb{Z}
$$

That is, if we add any integer multiple of $2\pi$ in the argument of a sinusoid, we get the same value out of the function (stated differently, sinusoids are $2\pi$ periodic). That means that we can drop the $2\pi kn$ out of the above expression for $x_2[n]$:

$$
x_2[n] = \cos\left(2\pi f_1 \frac{n}{f_s} + \phi_2\right)
$$

Now compare the two discrete-time signals again:

$$
x_1[n] = \cos\left(2\pi f_1 \frac{n}{f_s} + \phi_1\right) \\
x_2[n] = \cos\left(2\pi f_1 \frac{n}{f_s} + \phi_2\right)
$$

When we look at the sampled versions of $x_1(t)$ and $x_2(t)$, they appear to be sinusoids at the exact same frequency. In the discrete-time domain, since one or both of them were not sampled at at least the Nyquist rate, and the two frequencies have the relationship that was stated above, we are unable to discriminate between the two frequency components. This is what is called aliasing.


Cyclostationary signal intuition Images show a discussion I picked up from a PhD thesis about a cyclostationary process and need help interpreting it.

"In the time domain the upsampling process creates a signal whose distribution of amplitude is periodic or cyclic in nature" : I see what the upsampling has done but I fail to understand how the amplitude distribution is periodic?

The author goes on to say "yet within each of the  four  time  offsets  the  distribution  is  stationary,  leading  to  the  name cyclostationary" : is it stationary within the time offsets?

" In  the frequency domain four perfect images of the original spectrum are created due to the upsampling" how has upsampling caused this to happen?

I don't mean to discredit the authors work, this is merely a reflection of my lack of understanding.

I see what the upsampling has done but I fail to understand how the amplitude distribution is periodic?

The reason they can say this is because the symbols are assumed to be equally likely. Think about an upsampled BPSK stream, you'll have approximately the same amount of $+1$ pulses as $-1$ pulses so if you take a large enough block of samples it is periodic with period equal to the symbol period. In fact, this periodicity can be exploited to get an estimate of the symbol period (check out the paper Blind symbol rate estimation using autocorrelation and zero crossing detection). For enough symbols, the method does quite well for a blind estimator.

  The author goes on to say "yet within each of the four time offsets the distribution is stationary, leading to the name cyclostationary" : is it stationary within the time offsets?

Yes, as the author says: "within each of the four time offsets the distribution is stationary..." Cyclostationary means that the statistics change for different time offsets but has cycles. Say the peak of your pulse has a value equal to $1.0$. In noise, you sample the peak every symbol (cycle $=$ symbol period here) and you get a bunch of values: $1.1, 0.95, 1.04, 0.98, 0.89$ and so on. You can then collect statistics on these samples like the mean and variance. Cyclostationary means that changing the sample time by a half a symbol period will give you different statistics, but changing the sample time by a multiple of the symbol period (remember cyclo/cyclical) will give the same statistics. This is different from stationary process, where the time shift does not at all affect the statistics.

  "In the frequency domain four perfect images of the original spectrum are created due to the upsampling" how has upsampling caused this to happen?

You can see from the figure you posted, going from the plot labeled "Stationary Symbols" to "Upsampling" we are putting $3$ zeros between every symbol, or every symbol now occupies $4 \times$ as much samples (actual symbol + $3$ zeros = $4$ total samples). This is upsampling by a factor of $4$. When we do this operation, we get a compressed frequency domain axis and $4$ copies of the original spectrum (since we upsampled by $4$). I'm not deriving anything since this is a pretty standard DSP thing, check out
DSP Lecture 14: Continuous-time filtering with digital systems; upsampling and downsampling on Youtube starting at about the 56 minute mark for more in-depth upsampling explanation.


How do I plot the square of the amplitude response? I have calculated the transfer function of an FIR filter

$$
y[n] = x[n] + α · x[n − R]
$$
This is what I have
$$
H(z) = 1 + αz^{-R}
$$

Now I should plot the square of the amplitude response. So I calculated the frequency response by inserting e^jωTs for z. I should choose α = 1/2 and R = 3. So I get

$$
H(e^{jωTs}) = 1 + αe^{-RjωTs} = 1 +  \frac{1}{2} e^{-3jωTs}
$$

$$
|H(e^{jωTs})|^{2} = (1 + αe^{-RjωTs})^{²} = 1 +  \frac{1}{2} e^{-6jωTs}
$$

I am also not sure here about the algebra of the square of the magnitude.

Now I should plot from 0 to ωs and I think I am missing Ts. I should probably just choose a value, right? How is the usual approach?

If the impulse response $$h[n]=\delta[n]+\alpha \delta[n-R]\hspace{1cm}=\{1,0,0,...\alpha\}$$ with the ellipsis being R presumably, the Fourier transform $$H(\omega)=e^{-j 0 \omega}+\alpha e^{-jR \omega}=1+\alpha e^{-jR \omega}\hspace{1
cm}\exists\hspace{1
cm}\mathbb{C}$$
where magnitude is defined as $\sqrt{x^{2}+y^{2}}=|x+iy|$

Since we want the square of the amplitude response $=|H(\omega)|^2$, we can square the real and imaginary parts of the transform$$\Re{\{H(\omega)\}}^{2}+\Im{\{H(\omega)\}}^{2}\hspace{.5cm}=\hspace{.5cm}(1+\alpha \cos(-R\omega))^{2}+\alpha \sin(-R\omega)^2\hspace{.5cm}\exists\mathbb{R}$$
Then plot from $\omega=[0,\pi)$

Might look something like this with your parameters

Please tell me if this gives the expected response because I'm relatively new to continuous/algebraic Fourier theory


Weighted Nuclear Norm Minimization for Image Denoising Recently, I saw new published papers like

Shuhang Gu, Lei Zhang, Wangmeng Zuo, Xiangchu Feng, Weighted Nuclear Norm Minimization with Application to Image Demonising [pdf].

about denoising images using Weighted Nuclear Norm Minimization (WNNM) approach and I am wondering what the physical intuition behind it is.

The main idea is to subdivide the image into patches $Y_j$ then estimate the free-noise patch as the solution of

$$\min\limits_{X_j}{\frac{1}{\sigma_n^2}\|X_j-Y_j\|_F + \sum_i|w_i\sigma_i(X_j)|}$$

where $\sigma_n^2$ is the variance of the noise, $X_j$ is the $j^{th}$ denoised patch to estimate, $\sigma_i(X_j)$ is $i^{th}$ singular value of the matrix $X_j$ and the weights $w_i$'s are non-negative values and chosen in a non-ascending order to satisfy the convexity property.

More specifically, I would like to understand what that means if an image has the singular values of its patches (or regions) constrained to be sparse (according to the expression of the objective function) and I am wondering also if this optimization problem would yield characteristic patterns that are the same in all the images.

Most of the Denoisers in Image Processing make a simple assumption - The data has small number of freedom degrees while noise has high number.
Hence if we try to represent the given noisy data with small number of parameters we probably match the data while exclude most of the noise.

In the case above we try to limit the number of Singular Values of the base which represents the patches of the image.
The assumption means that most of the patches can be represented by a dictionary with small number of Singular Values -> Small number of degrees of freedom to create "Crazy Patches".


Conjugate symmetry of the DFT of real-valued sequences I have read about Fourier transformation that real signals are "mirrored" in the real and negative halves of the Fourier transform because of the nature of the Fourier transform.  For example, if we have the real signal $x$ whose length is $N$, its $X = FFT(x)$ should be conjugate symmetric. That means the $X^*[m] = X[N-m]$.
That is ok, and I have read it clearly HERE
My question: assuming you have a real signal $x = [-1, 1, -1, 1];$ its $X = FFT(x) = [0,     0,    -4 ,    0]$
So, in that case, we cannot say that $X^*[m] = X[N-m]$.  Why ?
Is what I said above is right ? in which cases is right?

We can! In your example, $N=4$, and the DFT is real-valued, so you get $X[k]=X[4-k]$, and that's true:
$$\begin{align}k=0:\;X[0]&=X[4]=0\\k=1:\;X[1]&=X[3]=0\\k=2:\;X[2]&=X[2]=-4\end{align}$$
Note that by definition the DFT is $N$-periodic, so $X[N]=X[0]$.


Polarity inversion in frequency domain (polar coordinates) I don't fully grasp how polarity inversion is performed in the polar coordinates of frequency domain. The frequency components of the signal do not change, so the amplitude part is going to stay whatever it is. That means the phase must change, but how? If I multiply the phase by -1 the time domain gets reversed. Phase shifting on the other hand only generates a delay...

The only thing I can think off is to reverse the amplitudes, but as far as I am aware the amplitudes are always positive in polar format. Am I wrong on this?

This is really basic, but here's how it works. If $X(\omega)$ is the Fourier transform of $x(t)$, then we obviously have

$$\mathcal{F}\{-x(t)\}=-\mathcal{F}\{x(t)\}=-X(\omega)$$

Furthermore,

$$X(\omega)=|X(\omega)|e^{j\phi(\omega)}$$

where $\phi(\omega)$ is the phase. From this we have

$$-X(\omega)=-|X(\omega)|e^{j\phi(\omega)}=|X(\omega)|e^{j(\phi(\omega)\pm\pi)}$$

because $e^{\pm j\pi}=-1$. So Jazzmaniac's comment that the phase needs to be shifed by $\pm \pi$ (or any odd multiple of $\pi$ for that matter) is of course correct.


What Is the Difference Between PCA and Karhunen Loeve (KL Transform)? I have been reading about Karhunen-Loeve or also known as KL transform and I see that when it is used to reduce dimension the procedure is identical to PCA, that is, for both methods the covariance matrix of the data is constructed and then the eigenvectors are calculated . I would like to know if there is any difference in this aspect that is being overlooked?

For discrete data both are the same - Finding set of orthogonal directions which maximizes the Variance (Energy) of data along them. Sometimes those are called the natural axis of the.
Since we're dealing with variance it is only natural both are calculated from the covariance matrix of data.

You may encounter places where KL might be even defined on the correlation matrix.

Anyhow, some nuances when dealing with them:

Usually when dealing with Continuous Function people use the term KL. PCA is preserved to discrete data (Though Mathematicians might work with the PCA in all cases as they have SVD for continuous data as well).
PCA is used for any discrete data where KL was specifically formed with random variables in mind.

In practice, when one has discrete data and tries to find the directions which maximizes the energy of the data (Dimensionality reduction, compression, etc...) both are the same.


Can I compute signal-to-noise ratio on a vector-valued or multi-channel signal? I have a dataset that contains multiple channels recorded over time (think EEG or EMG), and I'm working on a computational model for this data. When I have a working version of the model, I'd like to compare the model's predictions to some held-out test recordings. I thought it would be nice use some sort of signal-to-noise analysis to do this comparison, but I can't find much help on whether it even makes sense to try to compute SnR for multi-channel data.

Everything I've seen so far about computing SnR assumes that you have a single time-varying channel to analyze. Should I compute the SnR for each channel of my data independently, and then average the values across channels? How is this sort of analysis typically handled in signal processing?

You can take the way SNR is calculated in RGB images as a starting point.

In RGB images, the image SNR is calculated using the familiar formula

$$
\frac{\mu}{\sigma}
$$

Where $\mu$ is the "whole image" mean and $\sigma$ is the "whole image" standard deviation. To calculate these "whole image" values, each of the RGB channels is first converted to the common currency of luminosity by one of several formulas, generally weighted sums:
$$
L = W_{R}\mu_{R} + W_{G}\mu_{G} + W_{B}\mu_{G}
$$

Where $\mu$ are means and $W$ are respective weights for the R, G, and B channels of the image (in the Adobe formula, for example, weights are 0.2974, 0.6273, and 0.0753 respectively), and L is luminosity.

For whole image standard deviation, the variance includes the covariances of the channels, so the full formula is
$$
Var(L) = \sum_{i,j} (r_{i,j}^2Var_{R} + g_{i,j}^2Var_{G} + ...
$$
$$
b_{i,j}^2Var_{B}+ 2r_{i,j}g_{i,j}Cov_{R, G} + ...
$$
$$
2r_{i,j}b_{i,j}Cov_{R, B} + 2g_{i,j}b_{i,j}Cov_{G, B} )
$$
where $i$ and $j$ are pixel indices, which is a more expensive calculation and not as commonly used in industry. Oscar de Lama reports that the Adobe formula simply disregards the covariances:
$$
W_{R}^2Var_{R} + W_{G}^2Var_{G} + W_{B}^2Var_{B}
$$
and suggests a compromise formula:
$$
W_{R}Var_{R} + W_{G}Var_{G} + W_{B}Var_{B} ...
$$
$$
+ (2/3)(W_{R}\mu_{R}^2 + W_{G}\mu_{G}^2 + W_{B}\mu_{B}^2) ...
$$
$$
- (2/3)⋅(\hat{W}_{RG}⋅\mu{R}⋅\mu{G} + \hat{W}_{RG}⋅\mu{R}⋅\mu{B} + \hat{W}_{RB}⋅\mu{G}⋅\mu{B})
$$
Where $W$ are the earlier weights, and $\hat{W}(X,Y) = \sqrt{XY} / K$ where $K$ normalizes $\hat{W}$ such that $\sum \hat{W} = 1$.

If the vector-valued image weights all components equally, then W reduces to 1, simplifying the formulas.


Why Fourier series if Fourier transform can be calculated for both periodic and aperiodic? While learning about Fourier Transform after Fourier Series, That we can calculate Fourier transform of periodic signals too. If we can take the Fourier transform of periodic signal too then my question is why we Fourier series if Fourier transform can be calculated for both periodic and aperiodic?

You can't take the Fourier transform of a periodic signal, the integral diverges for all multiples of the period.

This can be handled by the theory of distributions, but the Fourier series is a better fit.


Increasing ADC output accuracy Say we have a 12-bit ADC with 30 MHz sampling frequency. I want to perform several DSP operations on output digital signal so that resulting signal would have increased accuracy to 14 bits and output sampling frequency would become 2 MHz. In nutshell I want a DSP block that converts a 12-bit ADC with 30 MHz Sampling rate to a 14-bit ADC with 2 MHz sampling. Even a slight tip for me to know where to look for would be helpful. (This project would be implemented on a FPGA.)

With typical sampling (not noise shaped) you get half a bit for every halving of the sampling rate (with proper filtering of the out of band noise prior to downsampling). So a decimate by two for example consists of a half band filter followed by a down-sample by two which is done simply by selecting every other sample.  Assuming white noise, the half band filter removed half the noise or 3 dB, and the the SNR of ADC’s goes as 6 dB/bit.
Therefore, if you went from 30 MHz to 7.5 MHz with a decimate by four, you would gain a full bit or 12 bits extended to 13 bits assuming noise is limited by quantization noise and the quantization noise process is sufficiently white as typically assumed and approximated.
It’s should be straight forward to see now how to go from 12 bits to 14 bits through successive decimations. (And the filtering and decimation can be done in one block each or distributed which gives flexibility in the filter designs used).
An additional caution when working with real hardware is to refer to the expected Effective Number of Bits (ENOB) for the actual ADC used and in what conditions of sampling, input frequency etc as the number of bits to be achieved and improved on, and pay attention to actual spur levels (as given by the Spurious Free Dynamic Range, SFDR) that will ultimately limit the SNR since decimation will do nothing to reduce the power level of in band spurs. Often these spurs are below the noise level due to quantization in which case dynamic range improvement can be gained through decimation.


Frequency shift property of Fourier transform Which one of the following is actually a definition of frequency shift property
$$e^{j\omega_0 t}x(t)\leftrightarrow X(j\omega - \omega_0) \tag{1}$$
$$e^{j\omega_0 t}x(t)\leftrightarrow X(j(\omega - \omega_0)) \tag{2}$$
Is frequency shift applicaple in all the situations?

Case I
$$1 \leftrightarrow 2\pi \delta(\omega)$$
if we apply frequency shift property we may obtain  $$e^{j\omega_0 t} \leftrightarrow 2\pi \delta(\omega -\omega_0)$$
which works according to result 2
Case II
$$u(t) \leftrightarrow \frac{1}{j\omega} +\pi \delta(\omega)$$
$$e^{-at} u(t)\leftrightarrow \frac{1}{a+j\omega}$$
which exactly isn't 1 or 2

Your Eq. $(2)$ is the frequency shift property; Eq. $(1)$ is wrong. This is easy to show:

$$\begin{align}\mathcal{F}\left\{e^{j\omega_0t}x(t)\right\}=&\int_{-\infty}^{\infty}x(t)e^{j\omega_0t}e^{-j\omega t}dt\\&=\int_{-\infty}^{\infty}x(t)e^{-j(\omega-\omega_0)t}dt\\&=X(j(\omega-\omega_0))
\end{align}$$

Your 'Case II' has nothing to do with the frequency shift property, because it is no frequency shift. For a frequency shift you need to multiply by a complex exponential $e^{j\omega_0t}$, and not by a real-valued exponential $e^{at}$.

You can come up with a rule for the Fourier transform of $x(t)e^{-at}$, $a\in\mathbb{R}$, but this is more tricky than the frequency shift property. For the frequency shift property, if you know that $X(j\omega)$ exists, then you know for sure that also the Fourier transform of $x(t)e^{j\omega_0t}$ exists. On the other hand, if $X(j\omega)$ exists, it is not certain that also the Fourier transform of $x(t)e^{-at}$ exists. This has everything to do with the Laplace transform and its region of convergence:

$$\begin{align}\mathcal{F}\left\{x(t)e^{-at}\right\}&=\int_{-\infty}^{\infty}x(t)e^{-at}e^{-j\omega t}dt\\&=\int_{-\infty}^{\infty}x(t)e^{-(a+j\omega)t}dt\stackrel{?}{=}X(a+j\omega)\end{align}$$

The last equality only holds if the integral converges, i.e. if $s=a$ is inside the region of convergence of the Laplace transform of $x(t)$, and if the Laplace transform of $x(t)$ has no singularities on the imaginary axis, i.e. if its Fourier transform has no Dirac delta impulses. The latter explains why the above 'rule' doesn't work for $x(t)=u(t)$ (as in your example).


DC gain from an impulse response I am studying for my exam in signal processing. For one of the old exam sets, a discrete impulse response of a filter is given as $h[n]$
\begin{bmatrix}
-3&0&-3&0
\end{bmatrix}
With a frequency response $H[k]$. What is the easiest way of calculating the DC gain from this information? I have tried using a DFT but did not see any usable conclusion.

The DC gain is simply the sum of filter taps or coefficients. This is the value of the frequency response at DC (i.e. $0\ \rm Hz$), or equivalently
$$
H(0) = \sum_{n = 0}^{N - 1}h[n]\tag{1}
$$
Because, for a digital FIR filter of length $N$ with impulse response given by Equation $(2)$
$$
\big\{h[n]\big\}, \quad\text{with}\quad 0\le n\le N -1\tag{2}
$$
the frequency response is the $z$-transform at the unit circle (i.e. $z = e^{j2\pi f}$), or in this case, equivalently
$$
H(e^{j\omega}) \equiv H(\omega) = \mathcal F \left\{h[n]\right\} \triangleq \sum_{n = 0}^{N - 1}h[n]e^{-j\omega n } = \sum_{n = 0}^{N - 1}h[n]e^{-j2\pi n f}\tag{3}
$$
Then you can think at DC (i.e. at $0\ \rm Hz$ or $f = 0$ or $\omega = 0$) to see the why of Equation $(1)$.
From DTFT to DFT
The DFT is the one practically computed in place of the DTFT. For finite-length signals (as is the case here), it provides the frequency-domain samples of the DTFT. Then $H[k]$, which is the DFT of the $h[n]$,  is computed for $N$-point as shown in Equation $(4)$
$$
H[k] = H(e^{j\omega}) \bigg\vert_{\omega = 2\pi k/N}\quad\text{with}\quad 0\leq k \leq N -  1\tag{4}
$$
i.e. $H[k]$ consists of equally-spaced (by $2\pi/N$) samples of $H(e^{j\omega})$.
Noting that
$$\omega \equiv 2\pi f \quad\text{with}\quad -\pi \leq \omega \leq \pi\quad\text{and}\quad -\frac 12 \leq f \leq \frac 12$$
With $\omega$ in [radians/sample] and $f$ in [cycles/sample].
In MATLAB
This can be done using MATLAB's freqz function as follows:
>> [h, w] = freqz([-3 0 -3 0], 1, [0 , pi/4])

h =

 -6.000000000000000 + 0.000000000000000i -3.000000000000000 + 3.000000000000000i

w =

                   0   0.785398163397448

>>

Giving you the flexibility to specify at whatever angular frequencies $\omega$  (in example above the evaluation is at two frequencies $0$ and $\pi/4$ rad/samples) you would want to evaluate the frequency response, (Here we're interested in the value at $\omega = 0$).  The magnitude response with these two points is shown below

Note that the two points in magenta correspond exactly to
>> 20*log10(abs(h))

ans =

  15.563025007672874  12.552725051033063

>>

The values at $\omega = 0$ and $\omega = \pi/4$ respectively.


How to organize MFCC features? I have computed 12d-MFCC features on 3 consecutive audio frames, say:

For 1st frame: a1 a2 a3 ... a12
For 2nd frame: b1 b2 b3 ... b12
For 3rd frame: c1 c2 c3 ... c12

How to organize these MFCC features, in a1 a2 ... a12 b1 b2 ... b12 c1 c2 ... c12 or a1 b1 c1 ... a12 b12 c12? Is there any difference between these two organization styles (for audio processing)?

Natural order would be to save it on frame by frame basis - that's the way you calculate them, otherwise you would need some gymnastics to do so. What's more, each 12 coefficients will be one vector in your 12D space - frame after frame.

Also it is easier to manage that in memory. In MATLAB like environments you simply use the matrices to store each vector in separate column. For C++, you can use array and array of pointers:

float **frames;
float  *data;
frames = new float*[num_frames];
data   = new float*[12 * num_frames];
memset(..)

After that you just set all frames pointers to first coefficient in each frame and you have easy way to access your features by using frames.

With interleaved features it wouldn't be so easy.


FFT Background subtraction in MATLAB Using MATLAB I need to subtract the background files from the real signal I am trying to measure. We took data of just the background noise prior to driving the system with some voltage, so we do not have to estimate the background, I need to know if i can subtract the two phase wise or how I can do this with the real and imaginary parts.

Thank you!!

I'm assuming since you didn't say that you're working with audio. But my answer generally applies if you're working with images/video as well. Unless the background you're trying to subtract is exactly the same, I doubt subtracting a phase-inverted version of the signal will work. It sounds like what you're trying to perform is spectral subtraction. This is a common technique where you simply subtract the average FFT magnitude spectrum from your signal. See this paper on spectral subtraction. If working with images/video, I assume that you could do this in 2 dimensions as well and the math would be similar.


How to determine if a filter is bandpass/stopband from its pole-zero diagram in z-domain How can we determine if a filter is bandpass or stopband, just by looking at its pole-zero diagram in z-domain?

For exmaple, if we have a system with third-order pole at the origin and a zero on the real axis at 0.5, what is the type of this filter?

In general, what is the effect of zeros outside the unit circle on frequency response magnitude?

For strongly selective filters that have well defined pass and stop bands and narrow transition bands and which, therefore, can approximate ideal brickwall filters to some desired degree, poles and zeros are generally distributed around and closer to the unit circle. Normally you would like to have all of them inside the unit circle but due to other reasons, such as linear phase FIR impulse response, some zeros can be put outside of unit circle too.

Effective frequency of a single pole/zero is the angle of its location from real axis counter clockwise. If that's a zero, it will stop frequencies nearby. If it's a pole it will amplify freqencies nearby.

By distributing those poles and zeros in groups you will be defining pass bands (that include band of poles) and stop bands (that include band of zeros), or bandpass filters as mixture.

When the magnitude of a single zero is away from unity, its selectivity decreases;i.e., pass-stop band distinction dissapears. When the magnitude is larger than one, it will also have an overall gain despite being a zero. If the magnitude is close to zero (zero at origin) then the gain approaches one.

Similarly, when the magnitude of a single pole is away from unity, its selectivity decreases;i.e., pass-stop band distinction dissapears. When the magnitude is larger than one, it will also have an overall attenuation despite being a pole. If the magnitude is close to zero (zero at origin) then the attenuation approaches one.


Minimum number of data points needed for a DFT to avoid spectrum leakage? Let $x(t) = 0.4 + 0.5 \cos (2 \pi f_1 t) + 2 \cos (2\pi f_2 t) + \sin(2 \pi f_2 t) + 1.5 \cos(2\pi f_3 t)$ where $f_1 = 3 kHz, f_2 = 5 kHz, f_3 = 8 kHz$. If $x[n]$ is obtained by sampling $x(t)$ with a sampling frequency of $f_s = 16 kHz$ what is the minimum number of data points needed for a DFT in order to resolve all the frequencies without spectrum leakage?

Can someone point me in the right direction on this?

UPDATE work:

With $f_s = 16000, T_s = \frac{1}{16000}$, we have:

\begin{align*}
  x(t) &= 0.4 + 0.5 \cos (2 \pi f_1 t) + 2 \cos (2\pi f_2 t) + \sin(2 \pi f_2 t) + 1.5 \cos(2\pi f_3 t) \\
  x[n] &= 0.4 + 0.5 \cos (2 \pi \frac{3}{16} n) + 2 \cos (2\pi \frac{5}{16} n) + \sin(2 \pi \frac{5}{16} n) + 1.5 \cos(2\pi \frac{1}{2} n) \\
\end{align*}

The discrete periods would be $\frac{16}{3}, \frac{16}{5}, 2$, right?

What you're actually asked to do is find the periods of the discrete-time sinusoids, and choose a DFT-length, such that all sinusoids occur an integer number of periods inside that window, hence no spectral leakage will occur.

I leave it as an exercise for you to determine the periods of the sampled sinusoids, but the end result of the DFT looks like this:

You can see that all frequencies are clearly resolved.

The problem (mentioned in the comments) of one sinusoid being located at exactly the Nyquist frequency (half the sampling rate) is not relevant for this toy example, even though in practice it is. What it means is that that sinusoid cannot generally be reconstructed uniquely. However, in your DFT you can see that there's a component at that frequency, and that's all they asked for.