Seven term Blackman Harris window I am currently trying to evaluate a high resolution ADC for which I need to window my output data for the FFT. Up until now I have been using the default 4 term Blackman harris window in MATLAB. However, I came across a few app notes that recommend using a seven term Blackman Harris for the SNR ranges that I am interested in.
I have tried looking quite a bit but I cannot seem to find an equation for implementing the seven term BH window, neither is it available in MATLAB. Can someone point me to a source where I can find this?
As far as I am aware Harris' classic paper on windows doesn't provide the equation either.

In this document you can find the coefficients of a seven-term Blackman Harris window. Ignoring the bizarre notation, it seems like the window is defined by

$$w[n]=\sum_{k=0}^6a_k\cos\left(\frac{2\pi kn}{N}\right),\quad n=0,1,\ldots,N-1\tag{1}$$

where $N$ is the window length and the coefficients $a_k$ are given by

$$\begin{align}
a_0&=0.27105140069342\\
a_1&= - 0.43329793923448\\a_2&=0.21812299954311\\a_3&= - 0.06592544638803\\a_4&=  0.01081174209837\\a_5&= - 0.00077658482522\\a_6&=  0.00001388721735\end{align}
$$

Unfortunately, they do not mention where those coefficients come from. Somebody must have run an optimization routine to minimize the sidelobe level for a given mainlobe width.

EDIT:

I think I found the relevant paper here. This is a corresponding Matlab function.


Is it possible to decompose an image into basis images using MATLAB? I have read that:

  Using Fourier decomposition any arbitrary image can be represented as
  summation of orthogonal basis images.

I want to see the basis images for any image say Lena or Cameraman, etc.

So can anybody help me to decompose an image into basis images and reconstruct back the same image from basis images using MATLAB?

It doesn't make any difference what image you choose. The basis images are the same. When decomposing a given image into basis images we find the coefficients only and define the given image as a weighted sum of all those basis images.

In the following, you can see the basis images based on 2D DFT of size $8$ ($64$ basis images in total)

Note the different 2-D frequencies. Any given $8\times8$ image is a weighted sum of these 64 images.


Waveform independent frequency analysis If a musician (or a vocal coach) hears a human voice at a certain note, he/she can easily find the same note on piano.
Not only for human voice, if the waveform is square, or created with a synth, they will be able to figure out which note it is.
However, using DFT (or FFT) gives mixed results for any waveform other than sinusoidal. It makes sense, because we use sinusoids as base functions (or base vectors).
Now, I want to learn if there is any way we can implement a signal processing algorithm in order to detect tones regardless of their waveform. Googling didn't help. Can you help me DSP SE?

What you are asking about is typically called "pitch detection". There are 100s of articles written on the topic. A good starter would be this one: https://ccrma.stanford.edu/~pdelac/154/m154paper.htm
If you need fairly high precision (for example a tuner application), I recommend a phase locked loop (PLL) or delay locked loop. There is a fundamental trade off between precision and tracking speed and with a PLL you can adjust this easily.


Problem getting DTW to work for speech recognition I used my ipad to record the words "This. Is. A. Test. This. Is. A. Test.", then downloaded it to my computer as a .MOV file and used ffmpeg to extract audio as WAV. I then used the tuneR package to load the WAV file into R.

library(tuneR)
sndObj <- readWave('myinput.wav')
plot(sndObj@left, type = "l")

It loads fine and the plot looks ok. The problem comes when I segment the words and calculate DTW distances between all pairs of words (using the DTW package in R).

I can't find any evidence that DTW distance is smaller between the correct pairs than between any other pairing (i.e "this->this" is about the same DTW distance as "this->a". Am I missing any important preprocessing?

Below: sample heat map of actual vs expected DTW distances.

Usually DTW includes mel-cepstrum feature extraction and cepstrum mean normalization. I'm not sure if you have it.


DFT of sum of sinusoids with random zeroed samples I have a noisy signal consists of sum of sinusoids, I had a situation where some of the time domain samples are zeroed randomly.
The following figure shows the DFT before and after zeroing, the number of zeroed samples are 40 out of 128 samples.

Obviously there is no much degradation of the DFT spectrum,.
What is the math behind this, and I would like to get some keywords to search any relative subject that discuss the limitations of zeroing.

The math is well known; it is the convolution theorem for the DFT.

In this specific case:

$$DFT\left\{f[n]\cdot z[n]\right\} = DFT\left\{f[n]\right\} * DFT\left\{z[n]\right\}$$

Where:

$f[n]$ denotes your original sequence of 128 samples

$z[n]$ denotes a sequence of 128 samples, each randomly selected from $\{0,1\}$

'$*$' denotes cyclic convolution, since we're talking about the DFT.

The results of this convolution are visible in your second plot, particularly at very positive and very negative frequency bins.


Building a FIR filter from an arbitrary frequency-magnitude response curve (eg. Least Squares fitting) I'm trying to create a FIR filter from a magnitude equation, where my starting equation provides the magnitude (amplitude) between 0 and 1 for any given frequency in Hz.

I posted the magnitude equations I would ideally like to convert into an FIR filter here:

https://math.stackexchange.com/questions/3666215/how-to-do-an-inverse-fourier-transform-of-this-equation

It seems there are more than one ways to get the result, but as I have never made a FIR filter before I'm not sure of the best way.

Option #1: Integral-Based Solution

One suggestion I was told was that I need to take the magnitude equation, do an Inverse Fourier Transform on it (which depends on the integrals), and then from that equation I can somehow construct a FIR filter.

But the integrals involve crazy operations like "Ei" and "erf" that probably aren't practical. So to do that, I would need to make a polynomial line of best fit for the magnitude equation and then use that for the reverse Fourier and FIR build.

Option #2: Least Squares Fitting

Alternatively, I read here a different suggestion for what I think is essentially the same problem:

Design a linear-phase FIR filter approximating the magnitude of a given IIR filter

That shares instructions on how to build a FIR filter to copy a IIR filter. If I understand, they used the "freqz" operation in Matlab to create a vector of frequency-vs-magnitude data points from the IIR filter. Those points were then fed into a least squares FIR filter design called "lslevin.m" which would then compute the FIR coefficients (or provide some output needed to get the coefficients).

If I'm reading that correctly, that approach would entirely solve my problem. I wouldn't even need the "freqz" operation, since I already have the amplitude/magnitudes at every frequency. I can just input those directly to the least squares FIR filter algorithm. It would also be linear phase which would be ideal.

Am I understanding this correctly, and if so, will that give me what I need?

Option #3: Directly Mapping Frequency-Amplitudes to FIR Bins

As a third option, I was recently told another approach:

  If you can get a value for the desired magnitude at any given
  frequency, the easiest approach is to sample the spectrum: simply
  compute an array of desired magnitudes with frequencies matching the
  FFT bin-frequencies for an FFT of your array size. No integration
  required. Essentially you just build the array that you would like to
  get out of FFT if you had the perfect filter.

  Then you can take IFFT to get a (periodic) signal that matches the
  desired spectrum at the compute bin-frequencies exactly. At this
  point, you want to apply some window function (with basically the same
  considerations as for any other windowed filter design), so you get a
  single IR (rather than a cycle of a periodic signal). This will smooth
  out the response slightly, so it won't be exact anymore, but assuming
  your original magnitude is reasonably smooth with respect to the
  length of the filter (ie. the frequency samples we took don't vary
  widely from one bin to the next) the windowed approximation is usually
  decent and you can make it better by computing a longer filter (by
  sampling more points and computing a longer IFFT).

  This type of design isn't "optimal" in any sense, but it's nice and
  simple and it works as long as you can tell the desired magnitude at
  any given frequency.

Is that going to be the same or different from the least squares method and in what ways? Which approach will be cheaper or more accurate?

Am I understanding this correctly, and if so, will that give me what I
  need?

Yes you are understanding everything correctly. If your "need" is a filter that will match the magnitude response then all choices should work. If you are interested in a very quick and move-on solution, then proceed with Frequency Sampling with a very large number of taps and move on. If you actually have a concern with efficiency (meaning you will be implementing the filter in hardware), then I recommend Least-Squares (although you could run into convergence issues). If your interest is education, then try all and compare the results yourself as nothing will beat the first-hand experience. All approaches will result in a solution.

The descriptions of the different approaches for designing FIR filters as stated by the OP are reasonably accurate, and given no constraint on total number of taps (and total delay) all could achieve the desired objective, so often the decision is just a trade of how much time you want to spend on designing the filter to achieve the best efficiency (least number of taps).

A fourth approach not mentioned is the minimax solution (equiripple) using the Remez exchange algorithm is also a common design approach.

The graphic below lays out my flow-chart for the typical steps in the FIR design process, with further details specific to estimating the order of the filter and the 4 major design approaches.

Estimating Order: Estimating the number of taps required in the FIR filter is determined before any further design efforts are made in any of the approaches. The number of taps required is driven from the derivative of the filter response versus frequency, both in magnitude and phase. The derivative of the magnitude response for typical low pass or high pass filter structures translates to very tight transition bands require a large number of taps, and fred harris' rule which is my favorite due to it's simplicity:

$$N \approx \frac{A}{22}\frac{2\pi}{\Delta\omega}$$

Where:
$A$: filter attenuation in dB
$\Delta\omega$: Transition band of filter in radians/sample ($2\pi = $ sampling rate).

The OP's case isn't a simple low pass or high pass filter, but this illuminates that the tighter the features are in frequency, the more taps will be needed.) Similarly the negative derivative of phase with respect to frequency is the group delay of the filter by definition. The OP does not care about phase but for general purposes we see how the slope of the phase vs frequency also drives the number of taps since the filter delay must be within the span of the filter.

The approaches are further detailed below by their common name, but as far as tap estimation, the approach above is reasonably close for the windowing, least squares, and equi-ripple design approaches. For frequency sampling a significantly higher number of taps may be needed  (could be 2x - 3x).

It would be common to start with a much higher number $N$ and then reduced the number of taps once the relative magnitude of the tails at the edge of the filter is reviewed. This can be reduced by pursuing the same algorithm with a smaller $N$, or truncating and windowing the solution achieved with a higher $N$.

Additional details on estimating the number of taps needed:

How many taps does an FIR filter need?

Design Coefficients:

The four common approaches for designing the FIR filter coefficients from target specifications are as follows:

Windowing (OP Option #1)
Frequency Sampling (OP Option #3)
Equiripple (Parks-McClellan, Remez Exchange Algorithm)
Least Squares (OP Option #2)

Windowing: (OP Option #1)

This method is to window the sampled continuous-time impulse response of the filter. The impulse response is determined from the inverse Fourier Transform of the frequency response, which is sampled and delayed to be causal and windowed to select the dominant response features based on the number of taps used.  Just selecting samples of the impulse response with no further modification is multiplying the impulse response by a rectangular window in time. Multiplication in time is convolution in frequency, so the desired frequency response will be convolved with a Sinc function in frequency. The Sinc function has the tightest main lobe for a given time span ($N$ samples), but has relatively high sidelobes. This results in the closest match of transitions but would have more ripple in the passband and stopband compared to specialized windows.

This is demonstrated below with an example of designing a perfect (target design, not achieved) low-pass filter. The impulse response of the low pass filter is a Sinc function, which is delayed in time in order to be causal (resulting in a linear negative phase slope versus frequency) and then truncated with a rectangular window to given the number of taps desired.

Once truncated, the other Sinc function in frequency due to the window (not the Sinc in time due to the desired filter response; in this case slightly confusing since two Sinc functions are involved) will convolve with the desired response, widening the bandwidth slightly and causing ripple in the passband and stopband. This is often improved by simply multiplying the desired impulse response with a different window shape that by design has very low side-lobes. The trade of doing this is a wider main-lobe which will serve to shift the exact locations of frequency transitions but will provide a much better match in areas where the frequency response is not changing rapidly. Increasing the number of taps resolves both cases.

Summary: Windowing is quick and easy if the inverse Fourier Transform of the target response is easily solved, and provides a reasonable result. It is sub-optimal in that other design approaches can achieve closer accuracy to target specifications with less number of taps.

Frequency Sampling: (OP Option #3)

This is the method I mostly discourage anyone from using (although I ultimately used that to demonstrate a quick approach for the OP's more challenging filter at his other post here: How to invert FFT these magnitude equations to get FIR filters?. showing that it does have utility if a quick and useful result is needed).

The benefit of this approach is that it is the simplest: You use the inverse DFT to create the coefficients of the filter. This will result in an exact match at the DFT sample locations, but the resulting response will have significantly higher ripple (deviation) for all frequencies in between.

A hybrid approach that would be very viable in this specific application of the OP's is to zero pad the high frequency portion of the target frequency response (which is the middle of the DFT array) allowing for a much longer duration for the impulse response (eliminating time domain aliasing). The resulting impulse response can then be windowed as in the window approach given above. Optionally the achieved impulse response can be resampled and then truncated/windowed. Basically we would be using a longer Inverse DFT to approximate the Inverse Fourier Transform when the result would be too difficult to compute directly (as in the OP's case).

When deriving the coefficients from the ifft, the frequency response needs to include the "negative frequency" components: for a filter with real coefficients this would be a mirror of the DFT samples from $0$ to $N/2-1$ at $N/2$ to $N-1$. Further, when the magnitude response is of only concern, the following will properly center the impulse response in the span of the filter as depicted in the diagrams above:

coeff = fftshift(ifft(ifftshift(mag_response)));

Here mag_response was generate with the frequency axis going from $-f_s/2$ to $+f_s/2$ so the inner ifftshift translates it to the format expected in the DFT going from $0$ to $f_s$, where $f_s$ is the sampling rate. (corresponding to $0$ to $N-1$). The final fftshift converts the noncausal result in the time domain (the coefficients are centered on $t=0$) to a causal filter by centering the coefficients in the center of the filter span.

Filtering by zeroing out FFT bins is the frequency sampling approach, which is further described here including why it is generally a bad idea (meaning will require a lot more taps than the optimal approaches equiripple and least-squares):

Why is it a bad idea to filter by zeroing out FFT bins?

Equiripple:

This is the go-to design algorithm when the desired accuracy result is constrained by peak limits. The equiripple design approach developed by Thomas Parks and James McClellan in 1972 https://en.wikipedia.org/wiki/Parks%E2%80%93McClellan_filter_design_algorithm uses the Remez-exchange algorithm resulting in an optimum filter design given a minimax constraint (minimize the maximum error). Optimum meaning the solution will be the mimimum number of taps given the peak error design constraint. It can be used to match arbitrary shapes and multi-band filters and is supported by functions in MATLAB (firpm), Octave (remez) and Python (scipy.signal.remez). For a more complicated filter it may have challenges converging which would be my only reason to pursue one of the earlier approaches above.

Least-Squares: (OP Option #2)

This is the go-to design algorithm when the desired accuracy result is constrained by root-mean-square (rms) limits. Given the performance of most applications where I need to use filters (wireless communications) is optimized by rms performance and not peak error, this is my filter design of choice. Like the equiripple design, this is an optimal filter design in that the solution will be the mimimum number of taps given the rms design constraint. Unless I was effectively given a "though shall not cross this line" constraint in the filter specifications, I would use least-squares.  It can be used to match arbitrary shapes and multi-band filters and is supported by functions in MATLAB/Octave (firls), and Python (scipy.signal.firls).  For a more complicated filter it may have challenges converging which would be my only reason to pursue one of the non-optimal design approaches above.

Further posts of interest:

FIR Filter Design: Window vs Parks McClellan and Least Squares


Variance of circular vector correlation Say I have two zero-mean vectors, with a size of N and the translation between them is k. Say the image signal std is $\sigma_s$ and the noise std is given by $\sigma_n$. What is the variance of the correlation value in case the correlation is estimated by Fourier transform, without zero-padding? Is there a way to find the variance for the Pearson-correlation too?

Presented are two vectors of length $l$ with a displacement of $k>0$ between them

$A_i=a_i+ω_i$

$B_i=b_i+ω_i$

where $ω_i\sim N(0,σ_n^2)$ and $a_i\sim N(0,σ_s^2)$. The displacement means that for $k<i<l-k$, $b_{i+k}=a_i$

Our estimation of the circular correlation, without zero-padding, takes the form of

$C_k=\frac{∑_{i=1}^NA_i B_{mod(i-k+m,m)}}{\sqrt{∑_{i=1}^NA_i^2∑_{i=1}^NB_i^2}} \frac{l}{l-k}$

Following derivation

$C_k=\frac{∑_{i=1}^N (a_i+ω_i)(b_{mod(i-k+m,m)}+ω_{mod(i-k+m,m)})}{\sqrt{∑_{i=1}^N(a_i+ω_i)^2∑_{i=1}^N(b_i+ω_i)^2}} \frac{l}{l-k}$

Expected value

The expected value for square of one sample is

$E[(a_i+ω_i )^2 ]=E[a_i^2+a_i ω_i+ω_i^2 ]=σ_s^2+σ_n^2$

Expected value for the circular correlation at the true peak location
$E[C_k]=\frac{σ_s^2}{σ_s^2+σ_n^2}$

which means that we use biased estimation of the true peak correlation.

Correlation value not at the true peak location is
$E[C_c]=0$

Variance

General expression for the variance of multiplication of two samples is given by
$var((a_i+ω_i)(b_i+ω_j))=var(a_i b_i+ω_i b_i+a_i ω_i+ω_i ω_j )$

For the correlated samples

$var(a_ib_i+ω_ib_i+a_iω_i+ω_iω_j)=var(2a_iω_i+ω_iω_j)=4σ_s^2 σ_n^2+σ_n^4$

For the uncorrelated samples

$var(a_ib_i+ω_ib_i+a_iω_i+ω_iω_j)=var(a_ib_i+2a_iω_i+ω_iω_j)=σ_s^4+4σ_s^2σ_n^2+σ_n^4$

Lets neglect the variance of the dynamometer. This leads to the variance of the true peak correlation value

$var(C_k )≈\left(k(σ_s^4+4σ_s^2 σ_n^2+σ_n^4 )+(l-k)(4σ_s^2 σ_n^2+σ_n^4)\right)\left(\frac{l}{l-k}\frac{1}{l(σ_s^2+σ_n^2)}\right)^2$

simplification leads to

$var(C_k)=\frac{kσ_s^4+l(4σ_s^2 σ_n^2+σ_n^4)}{(l-k)^2(σ_s^2+σ_n^2 )^2}$

and for a location which is not the true peak v≠k

$var(C_v)≈l(σ_s^4+4σ_s^2σ_n^2+σ_n^4)\left(\frac{l}{l-v}\frac{1}{l(σ_s^2+σ_n^2)}\right)^2$

here simplification leads to

$var(C_v)=\frac{l(σ_s^4+4σ_s^2σ_n^2+σ_n^4)}{(l-v)^2(σ_s^2+σ_n^2)^2}$

Therefore, the true peak value

$C_k\sim\left(\frac{σ_s^2}{σ_s^2+σ_n^2},\frac{kσ_s^4+l(4σ_s^2 σ_n^2+σ_n^4)}{(l-k)^2(σ_s^2+σ_n^2 )^2}\right)$

and the false peak value

$C_v\sim\left(0,\frac{l(σ_s^4+4σ_s^2σ_n^2+σ_n^4)}{(l-v)^2(σ_s^2+σ_n^2)^2}\right)$


Tuning a radio to a certain frequency I've a fundamental problem with understanding terminology used in DSP. What does it mean to tune a radio (for example, a software defined radio) to a certain carrier frequency $f$? The received signal is an arbitrary waveform and not a sine wave of frequency $f$. I understand AM modulation: you take a sine and manipulate its amplitude. The FM modulation I don't get: the modulated signal is not a sine, so why do we say that we tune a radio to a certain frequency in that case?

there are a lot of misconceptions in the question.

"The received signal is an arbitrary waveform and not [only] a sine wave of frequency $f$."

it's a collection of sine waves at different frequencies.  and "tuning" a radio or a filter or whatever is selecting frequency components in the neighborhood of $f$ while rejecting the others.


Frequency response of numerical derivative Analytical derivative of a function is equivalent to convolution of that function with $s$ in Laplace domain. Numerical derivatives are limited in bandwidth due to finite sampling rate, so they are not synonymous with convolving the signal with with $s$ term. At higher frequencies one would expect attenuation of the numerically differentiated signal from one that was computed analytically. Recently, I found that there are some differences at the low frequency limit as well which I cannot explain.

Attached is a plot of a signal sampled from a normal distribution (blue) and it's first derivative in time (red). As expected, at high frequencies the derivative signal begins to attenuate. But why does it not cross $\omega$ = 1 rad/s or 0.16 Hz as would be the case if the solution was obtained analytically?

Here's the code I am running in MATLAB

sr = 100000;
y = randn(1,sr);
dydt = y;
for i = 2:length(y)-1
    dydt(i) = (y(i+1)-y(i-1))*sr*2;
end
hold on, plot(abs(fft(y)));
plot(abs(fft(dydt)));
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')

When I try this, the results looks as expected. So you really have to explain in more detail what exactly it is that you're doing, because it doesn't seem to be a property of the discrete-time derivative.

EDIT: Now that I see your code, I'm convinced that the problem will disappear if you use

dydt = zeros(1,sr);

to initialize the derivative vector.


Why does fft require a power of 2 sample numbers? I'm a beginner in signal processing, I do not know if I'm right but what I understand so far is
if I have a signal of 44100 samples I have 22050 frequencies

so if I choose 512 samples for the fft the 43588 samples will not be analyzed
I can only edit the spectrum from 44100/512 = 86.1hz to 22050hz
if I want to edit the frequency 1hz I'll have to give 441000 samples to the fft but I can only give it or 65536 samples to analyze all the 22050 frequencies
then with 65536 given to the fft I can edit up to the frequency 32768hz

It doesn't.  An FFT is just a method or implementation of computing a DFT that runs faster if the length (N) can be factored.  It runs relatively faster (compared to a DFT of the same length) given a larger number of small prime factors of N, and is simpler for smaller numbers of unique prime factors of N (such as: only one).

Picking just one prime factor, such as the smallest, 2, makes textbook derivations and examples of FFT computations use less ink in the textbook and/or less chalk on the chalkboard.

But many modern FFT libraries are capable of computing an FFT for an N that is a composite product of several other (fairly small) prime factors, not just 2.


Best way to implement variable phase shift on FPGA? I have an incoming digitally sampled sine wave pulse, so my FPGA has the ADC level from an I and Q channel. I want to be able to shift the phase by some arbitrary amount.

What's the most efficient way to do this?

The most obvious way to shift phase would be just to multiply by $e^{i\phi}$, but that would be slow and wasteful. I've also seen someone recommend using an allpass filter with some phase shift.

Is there any other options I should be considering? Or is the allpass filter the best way to do this?

The allpass you're describing is a multiplication with multiple numbers. Hence, if you want to shift your phase, the only mathematically viable and easiest way is indeed a multiplication with $e^{i\phi}$.

I don't see why that would be overly slow – in a modern FPGA, you'd probably utilize DSP slices to do the (real-valued) two multiplications and additions that make up a complex $\cdot$ complex multiplication.

Now, if you need to increment your phase, so you need to *calculate the value of $e^{i\phi}$ first, then there might be room for optimization. In fact, digital tuning happens in FPGAs, in the DSP chains of software defined radio devices like the Ettus USRP series.

There, a pretty standard CORDIC algorithm was employed to approximate a complex sinusoid $e^{i2\pi f}$ of arbitrary relative frequency $f$. Have a look at the  USRP X3xx's GPL-licensed RX DSP DDC chain:

cordic_z24 #(.bitwidth(cwidth))
  cordic(.clock(clk), .reset(rst), .enable(run),
     .xi(to_cordic_i),. yi(to_cordic_q), .zi(phase[31:32-zwidth]),
    .xo(i_cordic),.yo(q_cordic),.zo() );

The cordic_z24 is really just a multi-stage CORDIC implementation with pre-computed atan coefficients.

It runs at complex 64MS/s even on 2005 Spartan FPGAs, and doesn't really blink much of an eye about 200MS/s on modern Kintex7 FPGAs.

It avoids explicit multiplication alltogether.


Dimension of channel matrix elements in MIMO systems In a MIMO system we have:
$y= \mathbf{H} \star \mathbf{x} +n$ with $x \in \mathbb{C}^{N_t \times 1}$ ($N_t$ the number of transmission antennas) and $\star$ the convolution product.

For the matrix, $\mathbf{H}$ we have: $\mathbf {H}=(h_{i,j}(t,\tau))$ with $i \in \{1,...,N_r \}$, $j \in \{1,...,N_t \}$, $t$ the time and $\tau$ the delay for the impulse response

Here is my question:
What is the dimension of $h_{i,j}(t,\tau)$?It's a number or a vector because in many book, there is $size(\mathbf{H})=N_r \times N_t$. But if $h_{i,j} (t,\tau)$ is a vector, it's false.

An example is article Survey of channel and radio propagation models for wireless MIMO systems.

Sorry, I'm starting again.
$\mathbf{H(t,\tau)}=(h_{ij}(t,\tau))$.

If I have $\mathbf{H}(t,\tau)=\mathbf{H(\tau)}$ and $N_t=1$ then $\mathbf{H}(\tau) \in \mathbb{C}^{N_r \times 1}$


Gaussian noise definitions I have an exercise that asks for us to determine the model for a system that goes like this:

  A battery manufacturer for cars is studying  how their batteries behave. The measurements are influenced by gaussian noise with center in 0V and standard deviation equal to 3 Volts. Nominal value for the Battery voltage is 14V. Determine the model for the system, as well as initial values to implement a Kalman Filter (in matlab).

It is quite clear what we want done. But I don't really understand how the centered values for the gaussian distributions could be different, and so I have no idea how to treat it. It is most likely that I have some concepts wrong, can someone explain how to deal with a noise that has a mean value different to the nominal value of the object in analysis?

The text:

  The measurements are influenced by gaussian noise with center in 0V and standard deviation equal to 3 Volts.

simply means that the Gaussian noise has $\mu = 0$ (the mean is zero) and $\sigma = 3$ (the standard deviation is 3 volts).

  Determine the model for the system, as well as initial values to implement a Kalman Filter

Without more information, a Kalman filter seems like overkill.

As @MBaz says, one model you could use is:
$$
r = V + n
$$
which has no time-dependent component. Making a Kalman-amenable model out of this might look like:
$$
V_{k+1} = V_k + m_k\tag{1}
$$
$$
r_k = V_k + n_k\tag{2}
$$
and we'd use the 14V nominal voltage to say $V_0 = 14$.

However, we are given no information about $m_k$ - how the true voltage might vary. We might also want to augment the state (perhaps the variation if $V_k$ depends on the derivative of $V+k$?) and so include a state transition matrix in front of $V_k$ in (1). If we did augment the state, we'd also probably want to add an output matrix on (2).

  Could you give me an insight about the initial covariance? Can it be reached from the information given?

Generally, I choose the initial state covariance to be VERY large. This is because I am generally given zero information about what it should be, as in this case. In this case, "VERY" might mean 30 (ten standard deviations).


Ridge finding in noisy 2D image I'm using MATLAB, and have experience with it but not with this kind of problem. I'm trying to identify white lines in noisy 2D images such as the one below. I can separate the brighter ones just by thresholding, but the dimmer ones (like in the bottom left) are tricky. Any suggestions? I guess I need some kind of ridge detection. It's worth noting that the background isn't entirely uniform... the brightness of the dim filaments might even be lower than the background noise in other parts of the image. So maybe some kind of local filtering instead of uniform background subtraction. As you can see I don't know where to begin.

Ideally I would will want to produce a skeleton describing all the lines for further analysis (measuring lengths). Of course this involves dealing with cross-overs and stuff too, but for now if I can just detect all the lines that would be a good start.

My proposal would be:

Filter noise
Create image without the lines
Subtract the image without the lines from the image with the lines
Apply threshold
Postprocess the image

Code:

close all;
% acquire input image
img = imread('ridges.bmp');
img = rgb2gray(img);
% remove some noise
img_gauss = imfilter(img,fspecial('gaussian',10,2));
figure();
imshow(img_gauss);
% remove all the lines from the image
img_filtered = medfilt2(img_gauss,[31,31],'symmetric');
figure();
imshow(img_filtered);
% subtract image without lines from image with lines
img_subtracted = uint8(double(img_gauss)-double(img_filtered)+127);
figure();
imshow(img_subtracted);
% threshold
img_thresh = img_subtracted>130;
figure();
imshow(img_thresh);
% use open to remove some noise
figure();
img_opened = imopen(img_thresh,ones(3));
imshow(img_opened);


Requirement of time step in convolution output? I am studying time domain analysis from book

signals and systems labotary with matlab

by alex palamides

when i was studying convolution command topic on page 187, author discussed use of time step with convolution command as shown in attached photo

Although comments and description written there but i couldn't understand its benefit/significance
Please kindly explain third rule especially
i am confused why he is using step there,why not go on without step?

The distinction is whether your x-axis is in units of time or units of samples, and if we want our digital system to remain constant with time independent of the sampling rate. Once our system is discrete we often prefer to work in (normalized) units of samples and not have to carry on the extra "step" term, in which case actual time will scale along with the sampling rate (this has been the more common approach in my experience, unless I am working on a mixed signal design).

A very easy way to see this is to consider a digital integrator (accumulator). We could follow the same process with convolution but using an integrator structure will be much easier, and clearer. First, the graphic below shows how a integrator which has the Laplace Transform $H(s) = \frac{1}{s}$ in the continuous time domain, can be mapped to $H(z) = \frac{T}{1-z^{-1}}$ in the discrete-time domain (using the mapping method of impulse invariance). Notice I stated specifically the discrete-time domain, meaning even though we are working with a discrete system, we still want to work in units of time (seconds). Therefore we must use the $T$ factor in the equation (which is "step" in your question).  If we work in units of samples, then T=1 and essentially disappears, simplifying our math.

So now the easy part: regardless of all the math and explanation above using "method of impulse invariance" etc let's just believe for a moment that the accumulator structure (keep adding the input to an accumulated value) is a discrete approximation of an integrator. If we add the factor T (or step) then the accumulator is a discrete-time approximation of an integrator.

Knowing that let's see it in action with reference to the graphic below showing the accumulator input and output.

If our horizontal axis is time (not samples) and if we want a system that keeps the same time performance independent of the sampling rate (meaning we achieve the exact same result with regards to time, the resulting waveforms throughout the system are just sampled faster or slower but the values vs time do not change), then in this case we must include the T factor.

Consider the simple case of the accumulator where the input is all ones. In discrete-time, regardless of the sampling rate an integration of a constant 1 at the input would grow to 1 after 1 second at the output. If the sampling rate was 10 Hz (for example) then T = 1/10Hz = 0.1. At 10 Hz, there would be 10 samples after 1 second, and without the T factor the output would grow to 10. By using T = 0.1, we get the result of 1 as expected at 1 second.  Change the sampling rate to 100 Hz, and after 1 second the output would grow to 100, but if we add the proper scaling T = 0.01, we still get the same result of 1 as expected at 1 second.

So in summary we see by including step (or "T" as in my figure) we produce results that are consistent with time (in seconds) regardless of the sampling rate. It is more common to work in units of samples (not seconds) and know that the system will scale along with the sampling rate, but using the units of time has convenience for mixed signal designs - or when a result must be reported in units of time.  We see this with normalized frequency for example where the frequency axis extends from 0 to 1 (cycles/sample) instead of 0 to $f_s$ (cycles/second) as described in this post:

What is normalized frequency


Understanding the Bilateral Filter - Neighbors and Sigma I have been reading about the Bilateral Filter in - Fast Bilateral Filtering
for the Display of High Dynamic- Range Images and have felt crippled and confused while understanding its working and the parameters.

In the equation for Bilateral Filter given as

How to determine the number of neighboring pixels?
Are all adjacent pixels the neighboring pixels? (top, right, bottom, left, diagonal)
What does it really mean when people say the following (underlined):

Please help me understand this. I feel lost as I try to understand and gain an intuition of the Bilateral Filter.

I plotted a graph of the Gaussian function (y-axis represents the value of Gaussian function;) with three-sigma values 0.2, 0.5 and 0.9 to help me understand but got nothing.

(blue = 0.2, green = 0.9)

I would start with the many resources on this site:

Is the Bilateral Filter a Solution of Some Variational Method?
How to Validate Bilateral Filter Implementation?
Comparison Between Guided Filter (Edge Preserving Filter) and Gaussian Filter.
What Is the Bilateral Filter Category: LPF, HPF, BPF or BSF?
Understanding the Parameters of the Bilateral Filter.

Regarding your questions, Let's address them one by one.

Q: How to Determine the Number of Neighboring Pixels?

The classic neighborhood used by Bilateral Filter is along the axis. So it is basically determined by Radius parameter. In order to make sense, the radius parameter is aligned with the $ {\sigma}_{s} $ parameter (Spatial Standard Deviation). Some examples of the connection between them can be found in:

How to Get Rid of Ripples from a Gradient Image of a Smoothed Image?
Parameters of Gaussian Kernel in the Context of Image Convolution.
How to Calculate Gaussian Kernel for a Small Support Size?
Gaussian Blur In MATLAB and Connection to Image Resolution.

Q: Are All Adjacent Pixels the Neighboring Pixels?

Well, Theoretically you can chose any way to set the neighborhood set (Circle, The whole image, etc...). But as I wrote above, the most common option it a rectangle (Square in most cases).

Q: What Are $ {\sigma}_{s} $ and $ {\sigma}_{r} $ For?

The parameters determines the range of values in the range and spatial domain.
Let's assume our image values are in the range [0, 1] and the spatial indices are integers. Assume $ {\sigma}_{r} = 0.1 $, then pixels with values which are farther than $ \sim 3 {\sigma}_{r} $ will have effective zero weight since $ {e}^{ -\frac{ {0.3}^{2} }{ 2 \cdot {0.1}^{2} } } = 0.0111 $ before normalization.

Same logic for the Spatial parameter $ {\sigma}_{s} $. If we set it to $ {\sigma}_{s} = 4 $ pixels which are more than 12 indices apart will have negligible weight. This is why we limit the radius of the neighborhood and match it to the spatial parameter as adding pixels which have zero weight is waste of computing resources.


What's the difference between the Hough and Radon transforms? I'm familiar with the Radon transform from learning about CT scans, but not the Hough transform.  Wikipedia says

  The (r,θ) plane is sometimes referred to as Hough space for the set of straight lines in two dimensions. This representation makes the Hough transform conceptually very close to the two-dimensional Radon transform. (They can be seen as different ways of looking at the same transform.[5])

Their output looks the same to me:

Wolfram Alpha: Radon
Wolfram Alpha: Hough

So I don't understand what the difference is.  Are they just the same thing seen in different ways?  What are the benefits of each different view?  Why aren't they combined into "the Hough-Radon transform"?

The Hough transform and the Radon transform are indeed very similar to each other and their relation can be loosely defined as the former being a discretized form of the latter.

The Radon transform is a mathematical integral transform, defined for continuous functions on $\mathbb{R}^n$ on hyperplanes in $\mathbb{R}^n$. The Hough transform, on the other hand, is inherently a discrete algorithm that detects lines (extendable to other shapes) in an image by polling and binning (or voting).

I think a reasonable analogy for the difference between the two would be like the difference between

calculating the characteristic function of a random variable as the Fourier transform of its probability density function (PDF) and
generating a random sequence, calculating its empirical PDF by histogram binning and then transforming it appropriately.

However, the Hough transform is a quick algorithm that can be prone to certain artifacts. Radon, being more mathematically sound, is more accurate but slower. You can in fact see the artifacts in your Hough transform example as vertical striations. Here's another quick example in Mathematica:

img = Import["http://i.stack.imgur.com/mODZj.gif"];
radon = Radon[img, Method -> "Radon"];
hough = Radon[img, Method -> "Hough"];
GraphicsRow[{#1, #2, ColorNegate@ImageDifference[#1, #2]} & @@ {radon,hough}]

The last image is really faint, even though I negated it to show the striations in dark color, but it is there. Tilting the monitor will help. You can click all figures for a larger image.

Part of the reason why the similarity between the two is not very well known is because different fields of science & engineering have historically used only one of these two for their needs. For example, in tomography (medical, seismic, etc.), microscopy, etc., Radon transform is perhaps used exclusively. I think the reason for this is that keeping artifacts to a minimum is of utmost importance (an artifact could be a misdiagnosed tumor). On the other hand, in image processing, computer vision, etc., it is the Hough transform that is used  because speed is primary.

You might find this article quite interesting and topical:

  M. van Ginkel, C. L. Luengo Hendriks and L. J. van Vliet, A short introduction to the Radon and Hough transforms and how they relate to each other, Quantitative Imaging Group, Imaging Science & Technology Department, TU Delft

The authors argue that although the two are very closely related (in their original definitions) and equivalent if you write the Hough transform as a continuous transform, the Radon has the advantage of being more intuitive and having a solid mathematical basis.

There is also the generalized Radon transform similar to the generalized Hough transform, which works with parametrized curves instead of lines. Here is a reference that deals with it:

  Toft, P. A., "Using the generalized Radon transform for detection of curves in noisy images", IEEE ICASSP-96, Vol. 4, 2219-2222 (1996)


Tied rank adjustment for Spearman Correlation I have two 16 sample vectors, a and b as listed below:

a = 345.86  349.83  257.84  438.43  195.03  0   341.39  0   0   0   0   0   0   0   0   0

b = 313.30  349.83  0   394.07  195.03  257.84  0   0   0   0   0   341.39  0   0   0   0

My aim is to calculate the Spearman correlation between the two vectors. I know that it can be done easily in MATLAB by using the function:

corr(a,b,'Type','Spearman');

However, my aim is to calculate this step by step, i.e.

Calculating the ranks of vectors a and  b
Applying the formula $C=1−\frac{6\Sigma(R_x−R_y)^2}{n(n^2−1)}$
Adjusting the tied ranks (most important part)

Firstly, using the function corr(a,b,'Type','Spearman') I get the correlation coefficient as 0.5845.

Now as to calculate the same step by step, I am using the tiedrank function to calculate the ranks for the vectors a and b. From The Mathworks documentation:

  [R,TIEADJ] = tiedrank(X) computes the ranks of the values in the vector X. If any X values are tied, tiedrank computes their average rank. The return value TIEADJ is an adjustment for ties required by the nonparametric tests signrank and ranksum, and for the computation of Spearman's rank correlation.

Using this for the above example, I get the following:

Rank_a = 14 15  12  16  11  5.5 13  5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5

Rank_b = 13 15  5.5 16  11  12  5.5 5.5 5.5 5.5 5.5 14  5.5 5.5 5.5 5.5

And the Tied Adjustments that I get:

TIEDADJ_a = 495

TIEDADJ_b = 495

So, later I calculate the differences of squares of the ranks between a and b i.e. Sum of (Rank_a(i) - Rank_b(i))^2 which is equal to 214.

I do not know how to take into account the TIEDADJ_a and TIEDADJ_b values in addition to the above summation of difference of ranks calculated, so that I can achieve the correlation coefficient calculated by the MATLAB function directly.

After finding the two rank vectors Rank_a and Rank_b, just calculate the correlation coefficient of ranks:

c = cov(Rank_a,Rank_b)/(std(Rank_a)*std(Rank_a))

From Wikipedia, you can only use
$$C=1−\frac{6\Sigma(R_x−R_y)^2}{n(n^2−1)}$$ if the ranks are distinct integer values (in this case are not).

  If ties are present in the data set, this equation yields incorrect results.


Estimation / Reconstruction of an Image from Its Missing Data 2D DFT Given the 2D DFT of an image i.e. a NxM matrix of complex numbers, with some missing lines (or even partial lines), considering we have zeros in the missing positions.

Any suggestions for an algorithm to "interpolate"/reconstruct the missing values?

For simplicity I will show an approach Id' use on 1D signal (A row of real world image).
You will be able to extend it and I will add few remarks on how you can even gain from having 2D data.

The general idea is as sketched in Estimate the Discrete Fourier Series of a Signal with Missing Samples. The trick here is to exploit prior information.
In our case we'll use the smoothness prior of images by utilizing the Total Variation prior.

Problem Formulation

So let't define our problem. We have the data to estimate $ x \in \mathbb{R}^{n} $ and the DFT samples $ y \in \mathbb{C}^{m} $ where $ n > m $. Let $ F \in \mathbb{C}^{m \times n} $ be the partial DFT matrix such that $ F x = y $, namely it transforms the to be estimated data into the given DFT samples.

This allows us to write the optimization problem as following:

$$ \arg \min_{x} \frac{1}{2} {\left\| F x - y \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$

Where $ D \in \mathbb{R}^{\left( n - 1 \right) \times n} $ is the Forward Finite Differences Operator (Numerical approximation of the Derivative).

By intuition, we're looking for a vector $ x $ which its partial DFT is similar to $ y $ and it obeys local smoothness as real world images do. The parameter $ \lambda $ balances between the two.

Solving this is pretty easy using Proximal Gradient Method or ADMM. But even straight forward Sub Gradient Descent Method will do the trick here with the Gradient given by:

$$ \frac{\partial}{\partial x} \frac{1}{2} {\left\| F x - y \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} = {F}^{H} \left( F x - y \right) + \lambda {D}^{T} \operatorname{sign} \left( D x \right) $$

Results

I used the Lena image and chose the 130th row:

I removed 3 random from the DFT of the row and estimated the row using that as a starting point.

The estimation after 1000 iterations:

The RMSE is about the value of 1 pixel Namely it is usually go unnoticed.

MATLAB Code

The MATLAB code is available at my StackExchange Signal Processing Q60119 GitHub Repository.

Pay attention one need to handle the fact the data to be estimated is Real while operations create Complex data. I used abs() in MATLAB to generate real data from the complex data.

How to Extend to 2D Data

Well, you could just do as above per row of the image. but since you have 2D information you should use that for better results.
The trick is to apply the Total Variation prior in both directions.

By the way, on noisy image, if you don't want to apply any smoothing solve this as Least Squares problem (Just do $ \hat{x} = {\left( {F}^{H} F \right)}^{-1} {F}^{H} y $).


Precision measurement of sine wave amplitude with ADC I want to measure amplitude of a sine wave input precisely with a limited resolution ADC.

As an example suppose that I have $1\textrm{ MHz}$ pure sine wave input to the $320\textrm{ Msps}$ $10$-bit ADC. I beleived there is many redundancy in data and with some signal processing I could get more Precision than $10$-bit.

Is there any way that I can do this without any change in circuit, like adding noise or other hardware change?

A sine wave has infinitesimally little bandwidth. By rotating, filtering appropriately and decimating, you can reduce the sample rate very much.

Each of these filtering operations is typically a summing operation, in which you "average" out noise (which isn't your main concern), but also get a more precise estimate for the amplitude.

Decimation in DSP is very commonly done. You'd probably want to do that anyway – 320 MS/s is really no fun to deal with, and you don't need that bandwidth.

Of course, you can also correlate with a synthesized sine, and measure the correlation coefficient to get the power/amplitude.

Other options are things like proper spectral estimators – there's a lot to choose from, including Welch's method, or Pisarenko-based approaches. Especially if your signal is noisy, these might be interesting, but it really depends a lot on what exactly you're measuring and how you model your noise.


What is the relationship between passivity of impulse response and passivity of the system? What I mean by passivity of impulse response is that the energy of impulse response for this system is less than the input impulse signal energy.

For example, in a discrete-time system, considering the impulse response, the input is Kronecker delta function $δ[n]$, and the output is the impulse response $h[n]$. Suppose that
$$
h[n]=\sqrt2/2 , n=1,2
$$
and $h[n]$ is 0 for any other components. In discrete-time system, We calculate the energy of the signal with
$$
Energy = \sum_{n=-\infty}^{\infty} |h(n)|^2
$$
Therefore, the energy of the impulse response is 1, and the input signal energy($δ[n]$) is also 1, which means the system is passive for this input signal.

Then considering another input signal $δ[n] + δ[n-1]$, the energy of this input is 2. The output can be derived by convolution of impulse response and this input signal. We get:
$$
y[n]=\sqrt2/2 , n=1,3
$$
$$
y[n]=\sqrt2 , n=2
$$
and $y[n]$ is 0 for any other components. The energy of this output is 3, larger than the input energy. Therefore, It seems that the passive impulse response doesn't necessarily  make a system passive. But if a system is passive, the impulse response of the system should always be passive. Am I right with this conclusion? or can anyone give a rigorous proof for this?

Let me clarify some definitions first. There is no such thing as a "passive impulse response", there are only passive systems. If you like, you can call the impulse response of a passive LTI system "passive", but that's not common usage. From what I understand from your question, you probably mean by a "passive impulse response" an impulse response $h[n]$ that satisfies

$$\sum_{n=-\infty}^{\infty}|h[n]|^2\le 1\tag{1}$$

i.e., an impulse response with energy not greater than $1$. However, as you've shown yourself in an example, this is not sufficient for the corresponding system to be passive. Let's use the following definition of passivity [1]:

  A discrete-time system is said to be passive if, for every finite energy input sequence $x[n]$, the output sequence $y[n]$ has, at most, the same energy, i.e.,
  $$\sum_{n=-\infty}^{\infty}|y[n]|^2\le\sum_{n=-\infty}^{\infty}|x[n]|^2\tag{2}$$

It is more straightforward to investigate the passivity of a linear and time-invariant discrete-time system in the frequency domain than in the time domain. Let $H(e^{j\omega})$ be the frequency response corresponding to the impulse response $h[n]$:

$$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}\tag{3}$$

Since the Fourier transforms of the input and output sequences are related by

$$Y(e^{j\omega})=H(e^{j\omega})X(e^{j\omega})\tag{4}$$

the condition

$$|Y(e^{j\omega})|^2\le |X(e^{j\omega})|^2\quad\text{for all }\omega\tag{5}$$

holds if and only if $H(e^{j\omega})$ satisfies

$$|H(e^{j\omega})|\le 1\quad\text{for all }\omega\tag{6}$$

Integrating $(5)$ gives

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}|Y(e^{j\omega})|^2\le \frac{1}{2\pi}\int_{-\pi}^{\pi}|X(e^{j\omega})|^2d\omega\tag{7}$$

from which it follows via Parseval's theorem that

$$\sum_{n=-\infty}^{\infty}|y[n]|^2\le\sum_{n=-\infty}^{\infty}|x[n]|^2\tag{8}$$

must hold. So it is not condition $(1)$ which makes a system passive, but condition $(6)$.

A causal and stable transfer function with real coefficients satisfying $(6)$ is called a bounded real transfer function [2]. This type of transfer function is important for the realization of low-sensitivity digital filters.

[1] S.K. Mitra, Digital Signal Processing, A Computer-Based Approach, 2nd ed., p. 70

[2] P.P. Vaidyanathan and Mitra, S.K., Low Passband Sensitivity Digital Filters: A Generalized Viewpoint and Synthesis Procedures, Proc. IEEE, vol. 72, no. 4, April 1984, pp. 404-423.


Performing DFT twice on an image. Why am I getting an inverted image? I was asked to perform DFT on an image twice as a part of my school assignment. Why am I getting a blurry inverted image when I perform DFT on an image twice? Sorry, I'm new to image processing and signal processing. Can someone help me explain this mathematically?

Here's the result I obtained.
PS: I used FFT algorithms to perform the DFT. (Python and numpy, np.fft.fft2)

[Edit]: I haven't shifted the DC value, F(0,0), to the center as usual, for the magnitude spectra given here. (I haven't used np.fft.fftshift)

Here's the exact code I used.

import numpy as np
from matplotlib import pyplot as plt

img = plt.imread("Q4.tif")
dft = np.fft.fft2(img)
dft_of_dft = np.fft.fft2(dft)

spectrum = np.abs(dft)
c = 255 / log(1 + np.amax(spectrum))
spectrum = c*np.log(spectrum)

double_dft_spectrum = np.abs(dft_of_dft)
c = 255 / log(1 + np.amax(double_dft_spectrum))
double_dft_spectrum = c*np.log(double_dft_spectrum)

p, axarr = plt.subplots(1, 3, figsize=(12,4))
axarr[0].imshow(img, cmap = 'gray')
axarr[0].set_title('Original Image')
axarr[1].imshow(spectrum, cmap = 'gray')
axarr[1].set_title('DFT')
axarr[2].imshow(double_dft_spectrum, cmap = 'gray')
axarr[2].set_title('Double DFT')
plt.show()

It's a DFT property that if you apply DFT twice to input data, you get the original signal flipped (circularly). Stated mathematically for 1D case:

$$ x[n] \xrightarrow{ N-DFT } X[k] $$
$$ X[k] \xrightarrow{ N-DFT } Y[k] = N x[-k] $$

similar result can be shown for 2D case. And as you can see, the resulting output is the flipped (horizontally and vertically) original data.

You can see mathematical derivation in other questions on this site. I remember having answered that before. Please check this site for applying DFT or FFT twice.

Such as this answer...


Laplace Transform of $f(t+a), a>0$ where $f(t)$ is not periodic For $a > 0$, is there any known representation of the Laplace transform of $f(t+a)$ in terms of the Laplace Transform of $f(t) $

Note: In my application, $f(t)$ is not a periodic function and the functional form of $f(t)$ is not actually known a-priori, because I have to couple it to another set of equations.

Let $s = \sigma + j\omega$, the inverse Laplace transform of $f(t+a)$ is given by
$$f(t+a) = \frac{1}{2\pi j} \int_{\sigma-j\infty}^{\sigma+j\infty} F(s)e^{s(t+a)} \mathrm{d}s = \frac{1}{2\pi j} \int_{\sigma-j\infty}^{\sigma+j\infty} F(s)e^{sa}e^{st} \mathrm{d}s.$$

Hence the bilateral Laplace transform of $f(t+a)$ is $F(s)e^{sa}$ where $F(s)$ is the Laplace transform of $f(t)$. For the unilateral case, see Matt L. answer.

This is sometimes called the shifting theorem. See Theorem 12.16 here.


Selectively improve the frequency resolution of FFT (DFT) I have a digitally sampled signal(sampling frequency = 4000 Hz) which contains analog frequencies from 0 to 1000 Hz.
I wish to have a frequency resolution of 100 Hz(like 0,100,200,300 Hz .. 1000 Hz) and I have to improve the frequency resolution at the 0 to 100 Hz range by 10 Hz. That is at the frequencies 0,10,20,30,... 100.
So the frequency components I wish to detect are [0,10,20,...,100,200,300,...,1000].
In other words, I want more frequency resolution at the 0-100Hz range and less resolution at 100-1000Hz range.
How this is possible through FFT or related methods.

You can compute any single FFT result bin using a complex Goertzel  algorithm or filter.

Realize however that if you want to compute a number of bins on the order of log N or greater (such as 10 extra bins with 10 Hz bandwidth and spacing), it may be faster to do the full longer FFT and throw away all the extra FFT result bins.

So just compute 2 FFTs of 2 different lengths and ignore all the unwanted results.


How to implement a high pass filter digitally to remove the dc offset from sensor data I'm trying to implement a high pass filter on sensor data to remove the dc offset. I'm aware I can subtract the mean of the data but I'm looking for a frequency domain approach. I've used the matlab's filterDesigner tool but I'm facing this problem:
You can see in the time domain output (last plot in blue) that the filter takes some time to settle. I want to minimize that and get similar results as mean subtraction (the plot in orange). Here are the filter specifications.

Can anyone help with me with what parameters to tweek to have no "settling time"? I'm going to implement this on esp32 so I need the order to be low.
I'm pretty new to dsp so please be kind.

The filter you desire is often called a “DC removal” or a “DC cancelation” filter. Unfortunately there’s no way to design a useful IIR DC removal filter that has no “settling time”. Your designed filter’s “transition region” is the lower passband frequency (50 Hz) minus the upper stopband frequency (10 Hz). So your designed filter’s transition region is 40 Hz.
If you design a “sharper” DC removal filter having a transition region of 20 Hz then that filter will have a longer settling time. If you design a DC removal filter having a wider transition region of 100 Hz (a lower “performance” filter) then that filter will have a shorter settling time than the filter you designed. So tryingengineer, you’re faced with a trade-off; the higher the performance of a filter (narrower transition region) the longer will be its settling time.
By the way, IIR filters have nonlinear phase responses. If you seek a computationally-efficient  DC removal filter that has a linear phase response, see the material at:
https://www.dsprelated.com/showarticle/58.php


CCDE processing vs. Scilab function I have used Scilab functions to produce a low-pass filter for an audio signal and the coefficients for the associated constant coefficient difference equation (CCDE).  I then produced filtered signals by running the Scilab filter() function and by running my implementation of the CCDE on the audio signal.  The results are identical.

The Scilab filter() function runs considerably faster, by perhaps a factor of 30.

I am new to DSP and I am trying to better understand what the Scilab filter() function is doing that allows it to use the same coefficients so efficiently.  In looking at the associated Scilab files, it looks like there is some compiled code behind the filter() function.

Any pointers as to technique or reference materials would be appreciated.

Scilab's filter, for short coefficient vectors, function implements a linear convolution in C code; that alone, since there's no python to actually be evaluated here, just multiplication and addition, is much much faster than writing something in a scripting language that can't 100% be just-in-time compiled.

For longer vectors, scilab implements fast convolution; ie. it exploits the fact that (circular) convolution in time domain corresponds to point-wise multiplication in (discrete) frequency domain, and uses zero-padding and saving of overlaps to emulate the linear convolution (which filtering represents) with that.

So, either way, use your libraries when doing signal processing! Aside from the convolution, there's other things that are generally faster if done via clever usage of library functionality: For example, whenever you have a loop that looks like

sum = 0
for a, b in zip(vectorA, vectorB):
    sum += a*b

you'd be far, far better of doing a dot product of the two vectors.

You have to consider this: Your CPU is very fast at doing basic math operations – often, it can do for example 8 multiply-and-accumulates (MAC) operations in a single step. Compared to that, parsing the structure of the (precompiled, even) for loop, building temporary python objects to hold the individual values for a and b, and overwriting the sum object, thus removing the old object and replacing it with a new one, leading to garbage collection and so on, is way way way way more work than just doing the maths. I like to put it like this:

  Imagine you're tasked with multiplying a lot of numbers between 0 and 10, but the numbers you need to multiply are written in text form in a book.
  Reading that book will take much, much longer than the multiplications

That's how it is to use dynamic languages to do basic math operations.


Normalization of an image in MATLAB If I have an image whose range is from $-18.3667\cdot10^{-5}$ to $9.3127$, how could I normalize it on the same gradient so that they fall between a specific range, say fro $0$ to $255$?

First, shift: put the minimum to $0$, by compensating the actual minimum $m=−18.3667⋅10^5$ for every pixel: $p\to p - m$. Now your pixels are between $0$ and a new maximum $M = 9.3127 - m$. Finally you want the final image in $[0,255]$. The second operation is scale: multiply by something, so that $0$ remains at $0$, and $M$ is cast to $255$. So you have to multiply by $255/M$.


Proof of $E[y(n)]=E[x(n)] \, \sum h(k)$ If $h(n)$ is the impulse response of the discrete LTI system, $x(n)$ is the white noise process with variance $\sigma_x^2$ and $y(n)$ is the output, prove that

  $$\mu_y = \mu_x \, \sum_{k=0}^\infty h(k)$$

  where $\mu$ represents the mean.

My approach:

$$\mu_y= \mathbb{E}[y(n)]
=\mathbb{E}[x(n) \star h(n)]$$
Where, '$ \star$' means convolution
Now how can I proceed?

The standard meaning of white noise includes an insistence (whether implicit or explicit) that the mean is $0$. Thus, what you want to prove is trivially true: since
$$Y[n] = \sum_{k=-\infty}^\infty h[n-k]X[k] = h\star X \big\vert_{n},$$ the linearity of expectation (the notion that $E[aX]=aE[X]$ and that the expectation of a sum is the sum of the expectations) gives us that
$$\mu_Y[n] = E\left[\sum_{k=-\infty}^\infty h[n-k]X[k]\right]
= \sum_{k=-\infty}^\infty h[n-k]E\big[X[k]\big]$$
where all the $E\big[X[k]\big]$ have value $0$, and so both sides have value $0$.

More generally, if the $X$ process has constant mean $\mu_X$, then
$$\mu_Y[n] = \sum_{k=-\infty}^\infty h[n-k]E\big[X[k]\big]
= \sum_{k=-\infty}^\infty h[n-k]\mu_X = \mu_X\sum_{k=-\infty}^\infty h[n-k]
= \mu_X \sum_{i=-\infty}^\infty h[i]$$ is also a constant.


Difference between spectrogram and a figure in an audio signal i wrote a function in MATLAB which will divide a audio waveform (a speech) into smaller windows with overlap .The output will be a matrix X, each column of which contains the windows which i export. Then i calculate the FFT of each column, i remove half of fft bins and i visualize the measure of each FFT frame as an image, using the "imagesc" command.Then i made a spectrogram of my signal and it gives me this figure. Im wondering if this two figures are the same or not. Ok it seems that maybe they are the same but why the colors for example are different? Or is there any other difference between this two figures?

The spectrogram is displayed on the log scale.


Whether $y(t) = x(t/3)$ is a memoryless system? I am a newbie in signal and system, and feeling confused about the memoryless property, it says at the definition, if a system gives output only depends on present input, then its memoryless. So if at $t=3$, $x(1)$ is needed, whether this count as a memoryless system? Or it only consider the time shift like $y(t) = x(t-1)$?

You are trying to mimic what I did at time $t$.  You'll be $y$ for now, say I'm $x$.

Now suppose that our are mimicking my actions  with your own modification pattern (called $f$), at the very same time. Then, your location, depending on mine, is:

$$y(t) = f(x(t))\,. $$

Here, you only need the current observation $x(t)$ (hence, no memory).

Now, suppose that you are trying to mimic what I do at time $t$, but only three times slower ($t/3$). What you will do at $t=3$ (or $y(3)$), I did it at $t = 3/3$, or at one second ($x(1)$). You need memory.


Geometric representation of a signal using basis functions I'm trying to solve the following problem where I need to express $x(t)$ in terms of the given $φ_i(t)$ functions. (It is proven that the $φ_i(t)$ functions are orthonormal).:

Here's what I tried:
In order to find the coefficients, the following definition of the geometric representation of a signal was used. (Simon Hagkin - Communication Systems)

Using equation 5.5 above
$$ x(t) = \sum_{j=1}^{3} x_{j} φ_j(t) $$
Referring to the equation 5.6 above, I have tried to find coefficients ($x_{j}$) of each basis function. However, all three coefficients became zero.
This is how it was solved for $x_{1} $:
$$ x_{1} =\int_{0}^{4} {x(t) φ_1(t)} dt $$
$$ x_{1}  =\int_{0}^{1} {-1 . (1/2)} dt +\int_{1}^{2} {1 . (1/2)} dt +  \int_{2}^{3} {1 . (-1/2)} dt + \int_{3}^{4} {-1 . (-1/2)} dt$$
$$ x_{1} = 0 $$
for $x_{2} $
$$ x_{2} =\int_{0}^{4} {x(t) φ_2(t)} dt $$
$$ x_{2}  =\int_{0}^{1} {-1 . (1/2)} dt +\int_{1}^{3} {1 . (1/2)} dt +  \int_{3}^{4} {-1 . (1/2)} dt $$
$$ x_{2} = 0 $$
for $x_{3} $
$$ x_{3} =\int_{0}^{4} {x(t) φ_3(t)} dt $$
$$ x_{3}  =\int_{0}^{1} {-1 . (1/2)} dt +\int_{1}^{2} {1 . (-1/2)} dt +  \int_{2}^{3} {1 . (1/2)} dt + \int_{3}^{4} {-1 . (-1/2)} dt$$
$$ x_{3} = 0 $$
What am I doing wrong?

Feels like a bit of a trick question. Your answer is indeed correct. $\Psi_1$ ...$\Psi_3$ are an orthonormal basis but it's an incomplete basis. Any signal you can construct with this basis has 4 degrees of freedom, but you only have three basis functions. In order to have a complete basis, you need 4 functions. Turns $x(t)$ (scaled properly) would make the 4th basis function, since it's orthogonal to the other three.


Angle of arrival (AOA) estimation using FMCW radar using MUSIC algorithm I am working with FMCW phased array radar with only upchirps. The range doppler matrix is obtained using the two dimensional fast Fourier transform on multiple chirps.
I want to obtain the angular spectrum per range velocity cell, so have taken range doppler bin for every antenna element and now I have a sample vector rdSnapShot (1 X Nant), where Nant is number of antenna elements present.

For this functionality I am using MATLAB's pmusic (MUSIC algorithm), my question is that as MATLAB's function pmusic requires correlation matrix, is it correct to just multiply the sample vector (1 X Nant) with its complex conjugate and obtain correlation matrix of (size Nant X Nant).

N.B. I have just one range doppler bin for each antenna element.

Providing the MATLAB code:

nTargets = 2;
Rxx_outerproduct = rdSnapShot*rdSnapShot';

[powerSpectrumMUSIC,idxFreq] =pmusic(Rxx_outerproduct ,nTargets,'corr');
% Adjusting the frequency for zero centered angle
idxFreq = idxFreq - max(idxFreq)./2;
% Shifting the response to obtain zero centeric response
powerSpectrumMUSIC=(powerSpectrumMUSIC);
% Converting the frequency scale to direction of arrival scale
idxAngle = asind(idxFreq/pi);
angResponseMUSIC = (db(powerSpectrumMUSIC));

Yes, it is proper/practical way to obtain correlation matrix. just multiply your sample with its conjugate to get auto-correlation and with conjugate of other antenna elements sample to get cross-correlation. arrange them in matrix to get correlation matrix($N$$X$$N$).
you can also look out for MVDR implementations in matlab which uses similar method for obtaining noise correlation matrix.


Bilinear transformation confusion Wikipedia says in bilinear transformation from \$s\$ domain to $z$ domain relation is

$$\boxed{s \longleftarrow \frac{2}{T}\frac{z-1}{z+1}}$$

But here this relation is given like this

$$\boxed{w=\frac{z+1}{z-1}}$$

What is difference between this two?

If you read the link which describes the second mapping
$$ w = \frac{z+1}{z-1}$$

  Thus inside of the unit circle in z-plane maps into the left half of w-plane and outside of the unit circle in z-plane maps into the right half of w-plane. Although w-plane seems to be similar to s-plane, quantitatively it is not same

It states that $w$ is not the $s$ plane.

Therefore, if your aim is in studying the bilinear transformation for digital filter design then you can safely use the first, well known, bilinear transformation relation, and its inverse:

$$\boxed{s \leftarrow \frac{2}{T}\frac{z-1}{z+1}}$$

Otherwise, the paper you linked is dealing with stability analysis of discrete-time systems using bilinear transform of its own selection...


whether the system is linear or not for the given problem Given the system:
$$y(t)=x(t+1)+x(t−1)$$
is the system linear?

For a system to be a linear first it should satisfy zero input and zero output. How can we calculate output at 0 input if the system depends on future or past or both? Please explain with steps.

How can we calculate output at 0 input if the system depends on future or past or both?

Well, surely zero input just means;
$$
x(t) = 0~~~~\forall t
$$
and the $\forall t$ means for all time: positive and negative.

Substituting that into the equation:
$$
y(t) = x(t+1) + x(t-1) = 0 + 0 = 0
$$

So the system is homogeneous.

Well, as @Dilip points out, this isn't sufficient for homogeneity: we need the output to be $a y(t)$ for all inputs $a x(t)$. The case above just looks at $a=0$.

The next question is: does it satisfy additivity?

If
$$
x_{\tt total}(t) = x_1(t) + x_2(t)
$$
then
$$
y_{\tt total}(t)  = x_{\tt total}(t+1) + x_{\tt total}(t-1)\\
= x_1(t+1) + x_2(t+1) + x_1(t-1) + x_2(t-1)\\
= y_1(t) + y_2(t)
$$
where
$$
y_1(t) = x_1(t+1) + x_1(t-1)\\
y_2(t) = x_2(t+1) + x_2(t-1)
$$
so it satisfies additivity also.


What Information Is Lost by Histogram Equalization? I cannot think of a way to restore an original image after histogram equalization. But this latter seems to "enhance" (= increase contrast of) the image's details. What information is then lost by histogram equalization, if any?

The expectation of information is called entropy. The loss of information can hence be understood as difference in entropy between source and processed image, assuming no random effect was added.

Together with my answer on why contrast is not an appropriate measure of entropy, this gives us the simple answer:

The loss in information is simply the number of bits necessary to count the different images that have the same processed image (i.e. after histogramm normalization).

Now, the math here gets pretty different whether you assume that pixels have real-valued intensities or discrete ones.

In the continuous case, you'd want to realize that histogram equalization applies a function to the pdf of the intensity, thereby adding negative differential entropy.

In the discrete case, multiple intensities of the source image will be mapped to a single intensity in the output: You can make statisitics of these ("empirical CDF") and from that derive how much information you lose.


Magnitude Spectrum of a Linear Chirp Signal I am trying to understand intuitively the fft of a signal that changes in frequency in time.
Suppose I have the fft of a linear chirp signal sampled at 1 kHz for 2 seconds. The instantaneous frequency is 100 Hz at t = 0 and crosses 200 Hz at t = 1 second.
Using matlab, I have the following magnitude spectrum plot:
t1 = 0;
t2 = 1;
Fs = 1; %1khZ
f0 = 100;
f1 = 200;
T = 2;

figure;
subplot(2,2,[1,2]);
t = 0:0.001:T;
y = chirp(t,f0,t2,f1); %generate chirp signal
Y=fft(y); %apply DFT on the signal
F=linspace(0,1e3,length(y)); %create x-axis for plotting the spectrum
plot(F,abs(Y)./length(y)*2); %plot the magnitude spectrum

I then truncate the signal by getting only the first 100 samples. Now I have:

I have some questions.

Are the plots correct?
Is it correct to say that by truncating, the ripples produced during transform are minimized hence the peak at 100Hz can be easily distinguished?

Are the plots correct?

Mostly. Two things to consider: since your input signal is real (as it not complex), the spectrum is conjugate symmetric. So typically you would only plot the first half from 0Hz to 500Hz. If you want to plot the whole thing, it's better to circulate it and plot it from -500Hz to +500Hz. In most cases you would also use a logarithmic scale (in dB) for the amplitude and not a linear scale.

Is it correct to say that by truncating, the ripples produced during transform are minimized hence the peak at 100Hz can be easily distinguished?

Not really. It looks like you truncated both the signal and the FFT length. At this point your frequency resolution is only 10 Hz, which means it's large compared to your sweep width. There are still plenty of ripples: you just can't see them because of your poor frequency resolution. There are only two discrete points in your "passband" and you can't make a lot of ripple with just two points.
If you don't like the frequency domain ripple, you can simply flatten it out in the frequency domain. That creates a little bit of amplitude ripple in the time domain, but it's ruler flat in the frequency domain. Which one you want to choose depends on what you want to do with the signal.

understand intuitively the fft of a signal that changes in frequency in time.

Intuition can easily fail you so it's a good idea to stay close to the math. The whole concept of instantaneous frequency only makes sense for signals that have very small bandwidth over a decent amount of time. For most signals that's not the case: while you can technically calculate it, the result is often meaningless. What's the instantaneous  frequency of the sum of two steady state sine waves with different frequencies ?


Hamming window for filtering unwanted terms I am trying to filter off two unwanted terms in this frequency spectrum (the spectrum on top).
First, I used a convolution between a rectangular window to make the unwanted terms equal to zero.

the first problem is when I tried to centralize the remaining part and multiplying it by a cos ( 2 * pi * f *t) I didn't get any result.

And my second question about hamming window is that in an article about a similar issue, a hamming and a blackman window were used. So, is it possible to use hamming window to filter off the unwanted terms? if yes, then how do I choose the starting and the final point of the hamming window to leave only the desired part! I really need help to understand how does it work and how to use it.

Below is some code (in scilab which is similar to but not the same as matlab) that attempts to do what you're saying.

The output is the following graph which has:

The original signal
The signal multiplied by the window to leave just the left-hand peak.
The signal modulated with a complex exponential to shift the remaining peak to the centre (DC).

// 25709
N = 1024;
dc = 1 + rand(1,N,'norm');
fc = 0.388729874;
a = 1+0.5*rand(1,N,'norm');
t = [0:N-1];
phi = 0.187238947;
x = dc + a.*sin(2*%pi*fc*t + phi);
X= fft(x);

f = fftshift((fftshift([0:N-1]) - N/2)/N);

clf;
subplot(311);
plot(f, fftshift(abs(X)))

subplot(312);
wind = [zeros(1,N-500) ones(1,500) ];
plot(f,fftshift(abs(X.*wind)));

subplot(313);
y = ifft(X.*wind);
yy = y.*exp(%i*2*%pi*fc*t);
plot(f,fftshift(abs(fft(yy))))


What is End Effect in signal processing? What does End Effect in signal processing basically means? I would like to see some reference book or papers if you suggest me.
Thanks

The end effect refers to where many filters and similar processes produce spurious results at the ends of the dataset.

There are several reasons why this can occur. The most common one is when using moving window based filters the filter is not well defined at the boundary. For example consider a filter with a window of size 3. At your first data point ($n=1$) this filter is not well defined as $n=0$ does not exist. In order to calculate a value at this point the most common approach is to pad your data (often with zeros) to extend the width. As this is not real data it will skew the results for these points.

I can't find much literature on it. The book Biosignal and Medical Image Processing has a little bit.


Time location of the DWT detail coefficients using MATLAB When performing the Discrete Wavelet Transform in MATLAB using the command DWT or WAVEDEC, what it the exact time or pseudo-time location of the DWT coefficients?

At each level the time series is decimated by 2, this means that the coefficients should be localized at $2^j \, T_s$, where $T_s$ is the sampling frequency. However, the application of a reconstruction filter of length $M$ makes the coefficient arrays longer. And each level a longer coefficient array is then fed into a new DWT function, which modifies further its length. The arrays are also properly extended at each level, according to a criterion specified by the user using the command DWTMODE.

I see that in some MATLAB examples, just the central part of the array is shown using the command WKEEP. Is this approach sufficiently accurate, meaning that the external coefficients (those that are not central) are really NOT significant?

If I do not discard some coefficients, I clearly see a huge delay if I assign them to $2^j \, T_s$. I believe that the external coefficients are however necessary for the reconstruction.

How many coefficients can be discarded at each level, if they can be?

This question has multiple facets  (after comments), so I will focus on the principal.

First, regarding coefficient localization: a discrete wavelet coefficient depends on several signal samples. The number of coefficients influenced by a single sample in a continuous wavelet representation typically depends on the properties of the mother wavelet and the signal regularity. This is illustrated on the following picture, with the modulus and the phase of a complex scalogram.

Coefficients at discontinuities (in the signal, the derivatives, etc.) spread in cones of influence. This is well-described in many wavelet books. The situation is even more complicated when you discretize the wavelet plane: one should project the samples with prefiltering, take care of the discrete wavelet symmetries and the level of redundancy. Remember for instance that the DWT is not shift-invariant. Hence, the mask could change a bit.

I thereby propose two methods:

one heuristic, based on the deterministic part of the data: build a simple template signal of what you want to detect (e.g. a discrete Dirac), perform your favorite discrete wavelet (redundant or not) over shifted versions, undo the shift scale-wise, combine the envelop of the scalograms and threshold them to keep the top values (as a percentage of the maximum amplitude). You can used it as a binary or weighted mask.
one more involved, based on the stochastic part of the data: it is possible to compute, or estimate, the second-order characteristics of a "random noise" (like a Gaussian distribution. The decay of the covariance matrix can serve to assess the influence of a noise sample in its neighborhood. There were many papers on that topic. We notably deployed this approach with our $M$-band dual tree wavelets: they are slightly redundant, and therefore there are correlations between scales and wavelet trees. This is described, as well as pointers to the relevant literature, in section III of Noise Covariance Properties in Dual-Tree Wavelet Decompositions.

The resulting "regions of influence" were later used in A Nonlinear Stein-Based Estimator for Multichannel Image Denoising: the shape of the mask (across scales and subbands) defines a Reference Observation Vector (ROV), on which we estimate the "denoised" coefficient, based on generalized thresholding expressions.

The above was used primarily for denoising, but similar reasoning could apply for adaptive filtering, restoration, segmentation, etc.


Python, resampling tools dont work as expected I've tried both using scikit.samplerate.resample and scipy.signal.resample. With neither of them can I achieve the result I want.

Scipy.samplerate:

file1 = '../machine_learning/voice_snippets_stepvoices/t10_9steps.wav'
file2 = '../random_data/t10_9steps.wav'

sr, y = scipy.io.wavfile.read(file1)

secs = len(y)/sr # Number of seconds in signal X

samps = secs*16000 # Number of samples to downsample
Y = scipy.signal.resample(y, samps)

scipy.io.wavfile.write(file2, 16000, Y)
#librosa.output.write_wav(file2, Y, 16000)

With both the outputs from scipy/librosa-write_wav() I can't listen to the file. So I assume the file is corrupted and this would be a result of the scipy.resample step. Yet I cannot work out why. I have no error logs to go off. I've tried plotting y and Y on wavplots and they are both the same except for different length x-axis(samples).

Any ideas?

Scikit.resample:

file = '../machine_learning/voice_snippets_stepvoices/t10_9steps.wav'
new_file = '../testing/random_data/t10_9steps.wav'
ratio = 16000/48000

sr, y = scipy.io.wavfile.read(file)
y_new = resample(y, ratio, 'sinc_best', True)
scipy.io.wavfile.write(new_file, 16000, y_new)

In this one I do get an output but it changes somewhat noisy audio of someone speaking to what sounds like Skrillex... (Joking. but look at the graphs).

I didn't expect the resampling to have this effect especially since scikit.resample is supposed to be the best resampling tool for python. (and yes it does apply anti-aliasing).

Anyone have any idea why these techniques aren't working?

Thanks a lot (:

EDIT:

Thanks @hulappa for the solution. It was related to the output I got from  scipy.io.wavfile.read. I used this code to solve:

 # scale to -1.0 -- 1.0 - https://stackoverflow.com/questions/2060628/reading-wav-files-in-python

if y.dtype == 'int16':
    nb_bits = 16 # -> 16-bit wav files
elif y.dtype == 'int32':
    nb_bits = 32 # -> 32-bit wav files
max_nb_bit = float(2 ** (nb_bits - 1))
samples = y / (max_nb_bit + 1.0) # samples is a numpy array of float representing

#--------------

Alas, this only worked for scikit.samplerate.resample(). I still get an unreadable file from using scipy.signal.resample().

I have tried running scipy.signal.resample() on a sine wave like so:

fs = 48000 # sample rate
f = 5 # the frequency of the signal
x = np.arange(fs) # the points on the x axis for plotting
# compute the value (amplitude) of the sin wave at the for each sample
y = [ np.sin(2*np.pi*f * (i/float(fs))) for i in x]

secs = len(y)/fs # Number of seconds in signal X
samps = secs*16000 # Number of samples to downsample
Y = scipy.signal.resample(y, samps)

file1 = '../random_data/sine48000.wav'
file2 = '../random_data/sine16000.wav'

scipy.io.wavfile.write(file2, 16000, Y)
scipy.io.wavfile.write(file1, 48000, np.array(y))

I cannot listen to both the resampled wav file or the sine wav file. I assume this is a careless error but I will put this here in case an answer helps future readers.

The problem is that if your wav file stores samples as integers (typically 16 bit signed), $\texttt{scipy.io.wavfile.read()}$ returns an array of integers. If this is the case, you want to convert your signal into a floating point representation before applying resampling. For a wav file which stores values as N bit signed integers, this can be done by:

// Convert y to float array, with values in the interval [-1.0,1.0]
y_float = float(y / (2.0**(N-1) + 1))


How to prove relationship between mutual information and differential entropy? I need to prove that mutual information given by

$$I(X;Y)=\int_{x,y}f(x,y) \log_2 \left( \frac{\left (f(x,y)\right)}{f(x) f(y)}\right) \, dx \, dy$$

is equivalent to $I(X;Y)=H(Y) - H(Y|X)$ I am proceeding like this

$$I(X;Y)=\int_{x,y} f(x) \frac{ f\left(x,y \right)}{f(x)} \log_2 (f(Y|X)) \, dx \, dy - \int_{x,y}f(x,y) \log_2(f(y)) \, dx \, dy$$

where first term gives me $H(X|Y)$ and second term gives me

$$\int_{x,y}f(y) f(x|y) \log_2(f(y)) \, dx \, dy $$

and I am taking

$$\int_{x} f\left (x|y \right) dx =1$$

I mean to say if this equation is right what does this physically symbolize that given $y$ there must be some $x$ that means given output there must be some input with probability $1$? Am I doing the derivation right ?

Your derivation is correct. Your second integral can also be written as

$$
\int_{x,y}f(x,y)\log_2(f(y))dxdy=\int_y\log_2(f(y))\left(\int_xf(x,y)dx\right) dy
$$

where the integral in parenthesis is equal to

$$
\int_xf(x,y)dx=\int_xf(x|y)f(y)dx=f(y)\int_xf(x|y)dx
$$

You can understand this integral in two ways:

Like your idea: $\int_xf(x|y)dx=1$, which means that regardless of what value $y$ has, $x$ must have some value, so the integral over all possibilities is 1. You can understand that also by seeing $f(x|y)$ as some distribution for $x$ (which is parametrized by $y$). And, the integral over the support of a distribution needs to be $1$.
$\int_xf(x,y)dx=f(y)$ is the marginalization of $f(x,y)$ regarding the variable $x$. With the integration, you remove the influence (i.e. marginalize out) of the variable $x$ from the probability density.


Does "keying on" a sine wave at a zero-crossing reduce its bandwidth? I understand that a pure sine wave of infinite duration occupies no bandwidth, i.e. it is only the modulation of a carrier that gives it sidebands. Does the exact timing of a sudden modulation make any difference?

For simple on-off keying of a signal, it looks like https://dsp.stackexchange.com/a/24604/18819 is giving an equation to determine the relative power away from a carrier:

  The spectrum of your windowed sinusoid is just a shifted sinc function

…and the magnitude of this function is based only on the length of the windowed signal. In all discussions I've seen (e.g. CW bandwidth described) they discuss things like keying envelopes and such to reduce the splattering of significant power away from the intended carrier frequency. They don't discuss the precise timing of the envelope.

I'm wondering: does it make any difference at which phase of its cycle I "key" a sine wave on or off?

The Wikipedia article on Zero-crossing circuitry explains:

  …the switching device will "wait" to switch on until the output AC wave reaches its next zero point. This is useful when sudden turn-on in the middle of a sine-wave half cycle could cause undesirable effects like high frequency spikes for which the circuit or the environment is not expected to handle gracefully

Is this just a practical concern (i.e. non-ideal components in a real-world circuit) or does the spectrum actually change if I make a "chopped-off sinewave" from 0º to 360º, vs. a section from 90º to 450º? What about from 0º to 180º vs. 90º to 270º?! [Update: to be clearer, maybe I'm really meaning in each case e.g. "from 90º to (360n+270)º" for some relatively large n so that the signal lasts for a bit between start/end.]

In both of those pairings, the signal which starts and ends at 0 magnitude certainly seems "cleaner" than a signal suddenly starting like in the last example at +1 and ending suddenly at -1. But can this difference be demonstrated mathematically? Each pair starts/ends the sine wave "instantaneously" and they both keep it running for the same duration, so are they equivalent as far as "splatter"?

It is instructive to have a look at the expression for the spectrum of a truncated sine wave to see how the choice of the window length affects the decay of the spectrum. Let's define a windowed sinusoid $s(t)$:

$$s(t)=\sin(2\pi f_0t)\cdot\text{rect}(t,T)=\begin{cases}\sin(2\pi f_0t),&-T/2<t<T/2\\0,& \text{otherwise}\end{cases}\tag{1}$$

The Fourier transform of $s(t)$ is the convolution of the transforms of the sinusoid and the rectangular window:

$$\begin{align}S(f)&=\frac{1}{2j}\left[\delta(f-f_0)-\delta(f+f_0)\right]\star T\text{sinc}(Tf)\\&=\frac{T}{2j}\left[\text{sinc}(T(f-f_0))-\text{sinc}(T(f+f_0))\right]\tag{2}\end{align}$$

With $\text{sinc}(x)=\sin(\pi x)/(\pi x)$ we get from $(2)$

$$\begin{align}S(t)&=\frac{T}{2 j}\left[\frac{\sin(\pi T(f-f_0))}{\pi T(f-f_0)}-\frac{\sin(\pi T(f+f_0))}{\pi T(f+f_0)}\right]\\&=\frac{1}{2\pi j}\frac{(f+f_0)\sin(\pi T(f-f_0))-(f-f_0)\sin(\pi T(f+f_0))}{f^2-f_0^2}\\=&\frac{1}{2\pi j(f^2-f_0^2)}\Big\{f\big[\sin(\pi T(f-f_0))-\sin(\pi T(f+f_0))\big]+f_0\big[\sin(\pi T(f-f_0))+\sin(\pi T(f+f_0))\big]\Big\}\tag{3}\end{align}$$

The spectrum $(3)$ decays as $\propto 1/f$ because we have $f^2$ in the denominator and $f$ in the first term of the numerator. If the window length $T$ is chosen such that an integer number of periods of the sinusoid are inside the window, then the limits of the window will be exactly at zero crossings of the sinusoid. This condition is equivalent to $Tf_0=n,\;n\in\mathbb{Z}$. With this condition the first term in the numerator of $(3)$ vanishes because the arguments of the two sinusoids differ by $2\pi n$, and, consequently, both sinusoids are identical. As a consequence, in that case the $f$ in the numerator disappears and the spectrum $S(f)$ will decay as $\propto 1/f^2$ instead of $\propto 1/f$.

In sum, if the windowing does not cause any discontinuities (i.e., the window limits occur at zero crossings of the sinusoid), the spectrum of the windowed sinusoid will decay as $\propto 1/f^2$. This is generally the case for continuous functions. Such functions can of course decay even faster, but only if they also have continuous derivatives, which is not the case for a windowed sinusoid. If the windowed sinusoid is discontinuous, i.e., the window limits do not occur at zero crossings, then the corresponding spectrum only decays as $\propto 1/f$.


Signal Processing in Amateur Radio Astronomy I am an amateur astronomer and I took recent interest in Radio astronomy and I am also new to Digital Signal Processing. My team is trying to build a radio telescope. I am interested in the data analysis part. One thing that is predominant in the signals that is received is the noise.

This is the data from the Radio Telescope of The University of Illinois recorded on July 4th 2011 and the telescope was pointing towards the Cygnus constellation:

http://radio.uindy.edu/data/2011/07041059-2.txt

Note: I am interested only in the third column of data

I made a plot of the data using R and that turns out to be something like this:

This is the plot of the Cygnus constellation inferred by the University:

I manually matched the above image with my plot and took out a portion of my plot:

It seemed similar to me. The chances are that my manual matching is wrong. How do I tackle with the situation ?

I have just started learning DSP and please suggest resources and pointers for DSP techniques so that I may apply in radio astronomy.

The signal processing operation that measures how "alike" two signals are is called cross-correlation.  For cross-correlation to work both signals must be sampled at the same rate.  If they aren't, you must resample one of them to the rate of the other.


Eye pattern construction and interpretation I have read some text about the eye pattern (or eye diagram), but I do not understand how should I read it.

The wikipedia definition of it is this one:

  In telecommunication, an eye pattern, also known as an eye diagram, is an oscilloscope display in which a digital signal from a receiver is repetitively sampled and applied to the vertical input, while the data rate is used to trigger the horizontal sweep. It is so called because, for several types of coding, the pattern looks like a series of eyes between a pair of rails. It is a tool for the evaluation of the combined effects of channel noise and intersymbol interference on the performance of a baseband pulse-transmission system. It is the synchronised superposition of all possible realisations of the signal of interest viewed within a particular signaling interval.

So my first question is: how does an ideal eye pattern look like?

My second question is about how should I read it, and I think that in order to understand this I should know how it is built. For instance this page shows it in this way:

According to wikipedia, it is sufficient to make the superposition of all the signals represented in picture. I understand that the resulting signal is that shown in the last graph. But, how should I read it? Which are the advantages of analyzing it instead of the single signals properties (edges etc)?

It is the objective of the receiver to make the best estimate for each symbol as to what was transmitted. This is often done by ultimately determining a decision time in each sample (through timing recovery) on the waveform after it has been processed by the receiver (equalization and matched filtering) in which to sample the waveform and make a decision as to what symbol was transmitted.

The eye diagram superimposes the waveform for each symbol at the point in the receiver just prior to decision. The eye diagram readily shows the ability to make a decision between signal levels; so the ideal eye would show that there is a lot of margin both in the vertical and horizontal axis to allow for a minimum error rate (what we would call a wide eye opening). As the opening narrows, it would show that there is an increased sensitivity to noise since any added noise would cause a value to pass a decision threshold resulting in an error (this could be phase noise or jitter in the vertical direction and amplitude noise in the horizontal direction, or intersymbol interference which effects both directions).

I have this graphic pasted below that I created for one of my courses that may be useful in answering your question. This particular eye diagram is over two symbol durations. The waveform directly below the eye-diagram plot is replicated in the eye diagram as the darker blue lines, amidst the rest of the trajectories of the waveform from symbol to symbol over a much longer duration as given in the lighter blue lines:

Other examples of eye diagrams are shown below. The first is for QPSK, where we see the constellation on the left with the ideal sample locations in red, and the remaining trajectories from symbol to symbol everywhere else. On the right is the eye diagram formed from the real and imaginary components of the waveform.

And below is an eye diagram for a 16 QAM signal before and after the root-raised-cosine filter in the receiver. Here we clearly see the utility of the eye-diagram in that we see due to the proper filtering we have significantly increased our margin against noise. This would be similar to using the eye diagram to assess the effectiveness of equalization or to determine if certain filtering in our receiver is degrading our ability to obtain a low error rate due to the eye closing.


Zero Phase - Moving Average Filter How can I obtain a Zero Phase Moving Average Filter? I read for example in Matlab that filtfilt give you zero phase doing forward and backward filter, I dont understand well how that work, I think taking the same number of values in the past and future can give me that, but it isnt causal right?

A zero phase moving average filter is an odd length FIR filter with coefficients

$$h[n]=\begin{cases}\frac{1}{N},&\quad -(N-1)/2\le n\le (N-1)/2\\0,&\quad\text{otherwise}\end{cases}$$

where $N$ is the (odd) filter length. Since $h[n]$ has non-zero values for $n<0$, it is not causal, and consequently, it can only be implemented by adding a delay, i.e. by making it causal.

Note that you can't simply use Matlab's filtfilt function with that filter because even though you would get zero phase (with a delay), the magnitude of the filter's transfer function gets squared, corresponding to a triangular impulse response (i.e., input samples further away from the current sample receive less weight).

This answer explains in more detail what filtfilt does.

See also this related question.


Why zig zag manner scan is used in DCT for image compression? I am not sure why zig zag manner scan is used in DCT. I can guess that might because it moves from low frequency component to high frequency component and as most of the energy is stored in low frequency component that's why this way of scan is favored but not quite sure.

You are right! When it comes to quantization, bit allocation or issues like that scanning in zig zag manner seems to be most efficient.

As you said, (and I emphasize) in Natural Signals, lower frequency coefficients contain most of the energy. In a DCT matrix, as you go along a row the frequency is increasing in the first dimension and as you go along a column the frequency is increasing in the second dimension. A high frequency in any dimension probably contains less energy. So for example the
$(2,2)$ entry in an $8\times 8$ DCT matrix $\bf{usually}$ contains more energy than the entry $(1,8)$ in the same matrix.


Autocorrelation - Understanding reduced correlation at periodic time shifts using np.correlate (versus statistical autocorrelation)? I'm going through the Think DSP by Allen B. Downey, and I'm struggling to understand a specific aspect of np.correlate and how it differs from statistic autocorrelation. The question is at the bottom, but to give it here: why does the magnitude of autocorrelation drop off at periodic intervals of time shift using np.correlate but not the statistical autocorrelation function (pearson product)?

The presented signal is

Using statistical autocorrelation, which I believe normalizes the signal and unbiases it, it looks like:

Using np.correlate(signal, signal, "same") I get

The book mentions that "you’ll notice that the correlations computed by np.correlate get smaller as the lags increase. That’s because np.correlate uses the unstandardized definition of correlation; as the lag gets bigger, the number of points in the overlap between the two signals gets smaller, so the magnitude of the correlations decreases." - This is the part I do not understand. It's effectively a periodic signal, so why does the number of overlapping points get smaller with greater time shift? At some point the shift matches the period, and the correlation should be near 30 again? Does it have something to do with the lack of normalization in the latter?

The difference here is that statistical autocorrelation assumes a stationary power signal (so basically an infinite periodic signal) and does a normalisation to $[-1,1]$ and that autocorrelation just takes the finite signal as it is and does the classical convolution of the signal with itself, leading to just one point of perfect overlap and decreasing amplitudes before and after.