Bandwidth with complex sampling On the transmit side, I have a 20 MHz carrier frequency carrying a signal with bandwidth of 40 MHz (so 0 Hz to 40 MHz, center at 20 MHz).

On the receive side, I have a dual channel ADC with each channel sampling at 40 MS/s (channel 1 samples at 40 MS/s, channel 2 samples at 40 MS/s). The hardware is setup to input I and Q into the two channels.

Am I able to properly sample the received signal, without aliasing? I saw that this was discussed before ("Complex sampling" can break Nyquist?), and the answer seemed to be "BW = fs because you get the negative frequencies", and I'm not 100% sure if the "negative frequencies" thing is a caveat, or I can just do the math as if I am getting 80 MS/s (which I technically am), so I drummed up this numerical example.

Complex sampling does not "break" Nyquist.  IQ quadrature sampling produces twice as many bits per second of information (at the same sample rate for real or complex samples), and the 90 degree phase offset between the I and Q channel in those bits provides extra information about the spectrum.
One typical example to demonstrate aliasing is that the samples of DC and the samples of a sinusoid of a frequency at the sample rate look the same.  A single channel of samples at 40 MHz would alias 0 Hz and a 40 MHz sinusoid together.  Only if you strictly bandlimit the input spectrum to below Fs/2 in bandwidth would you prevent this aliasing.  If you feed a strictly real data vector to an DFT, the result would be conjugate symmetric, half the FFT result would be redundant, and the upper half can be called the negative frequency image.
But with quadrature sampling (using a 40 MHz sampling clock plus a 90 degree phased shifted copy), a DC signal and a 40 MHz sinusoid would be clearly differentiated.  I and Q would be the same for DC, but different for 40 MHz, e.g. if I is the same, the Q would be different, or vice versa.  Therefore the spectrum between half the IQ sample rate and the sample rate will not be aliased with spectrum below, as is the case with single channel sampling.  If you feed a DFT using the IQ data as a complex input vector, the lower and upper half of the result can be non-symmetric, thus represent independant information, thus twice the bandwidth (between 0 and Fs), compared to a strictly real DFT of data at the same sample rate.


Implementing discrete Poisson equation wtih Neumann boundary condition I am trying to reconstruct an image from gradients in an arbitrary shaped region of an image.  I understand how to implement a discrete 2D poisson solution with Dirchlet boundary conditions. Using http://en.wikipedia.org/wiki/Discrete_Poisson_equation#On_a_two-dimensional_rectangular_grid , you just replace any of the u_i,j that are on the source side of the boundary with the boundary values that you have. However, I don't understand how to implement the Neumann boundary condition. If you want to "match" the external derivative, wouldn't you have to just set the nodes that are on the target side of the boundary to a value that makes the derivative between the source boundary pixel and the target boundary pixel equal the derivative outside the boundary? If you do this, haven't you just simply "filled in" one node around the boundary and provided a slightly smaller version of the Dirchlet problem?

Can anyone explain this?

For image in painting and blending a pixel its new value is dependent on its neighboring pixels. These neighbors may or may not lie on a border. The aim is to interpolate and create smooth transitions. It is thus not only dependent on a border condition, nor does a value on the border get replicated exactly to inside pixels.  Should the area between boundaries be exactly 1 pixel wide, this pixel will most likely be the average of the boundary values.

A simple blender algorithm will work like this. Given:

A background image $\textbf{b}$
A foreground image $\textbf{f}$
A mask $\textbf{m}$

where we will put the foreground image $\textbf{f}$ onto the background $\textbf{g}$, where the mask $\textbf{m}$ is true.

To avoid a visible seam we will make a correction image $\textbf{c}$. The run-of-the-mill patching to form new image $\textbf{g}$ then is:

$g(x, y) = \left \{\begin{array}{cc}
  b(x, y) & \textrm{if } m(x, y) = 0 \\
  f(x, y) + c(x, y) & \textrm{if } m(x, y) = 1
 \end{array} \right .$

This image $\textbf{c}$ is obviously some form of interpolation, but informally we want two things from $\textbf{c}$:

$\textbf{f}$+$\textbf{c}$ and $\textbf{b}$ should match on the edges
It should be 'smooth' as to not create disturbances

Smooth transitions can be attained by simply stipulating that a pixel its new value is the weighted average of its four neighbors (up, down, left and right):

$c(x,y) = \frac{c(x-1, y)+c(x+1, y)+c(x, y-1)+c(x-1, y+1)}{4}$

And for boundary pixels the corresponding values of $\textbf{c}$ are replaced by $\textbf{b}-\textbf{f}$. If for instance for $c(x,y)$ the right neighbour is on de boundary:

$c(x,y) = \frac{c(x-1, y)+c(x+1, y)+c(x, y-1)+(b(x-1, y+1)-f(x-1, y+1))}{4}$

This means that $\textbf{every}$ pixel of $\textbf{f}$ under mask $\textbf{m}$ gives a linear equation, and each pixel in the patch is also a variable.

The result is a square linear system:  $\textbf{A}\textbf{c} = \textbf{d}$, which can be solved for  $\textbf{c}$.

Note that

For boundary pixels, the values of $\textbf{b}$ and $\textbf{f}$ are $\textbf{constants}$ in this system, not variables. They thus go `to the right' of the equals sign.
Pixels outside the masked region, as defined by $\textbf{m}$ are not variables.

See also

J. Jia, J. Sun, C.K. Tang, and H.Y. Shum. Drag-and-drop pasting. In ACM
Transactions on Graphics (TOG), volume 25, pages 631-637. ACM, 2006.
P. Perez, M. Gangnet, and A. Blake. Poisson image editing. In ACM
Transactions on Graphics (TOG), volume 22, pages 313-318. ACM, 2003.


Correction of signal after High-pass RC-filter I have one question about digital correction of signal after RC filter which is high-pass.

Let me explain detail.

I have one simple signal conditioner. It has RC filter on the signal input with cut-off frequency equals 1.3 Hz. It means that after this RC filter my signal has lower amplitude from 1.3 Hz to approximately 15 Hz. After ADC I get this signal on my computer in digital form. After that I use Fourier transform.

Unfortunately, I need correct amplitude information in range 2-20Hz. So I need correct amplitude of signal in that range by digital transformations.

Can I do this? Is it normal practice?

Yes, this seems like a normal thing to do. However, note that the signal-to-noise ratio will be degraded in the 1.3 Hz to 15 Hz band. Your signal will be smaller but the quantization noise of the ADC will stay the same.


DFT of time reversed signal I was looking into proof and find something strange:

The last part we obtain from DFT definition.
$$X[k] = \sum^{N-1}_{n=0}x[n]W^{kn}_N, \quad\text{Where}\quad W^{kn}_N = e^{-j\frac{2\pi}{N}nk}$$

So, looking at $W^{-k}_NX[N-k]$ we can see that when $k = 0$ we will have $W^{0}_NX[N]$ and we know that our DFT coefficients only can be from $[0,N-1]$.

So my question is: what does it mean when we have $X[N]$ DFT coefficient? Or did i understand something wrong?

You can evaluate $X[k]$ (and $x[n]$) outside the domain $[0,N-1]$. The sequences are periodic with period $N$, so $X[N]=X[0]$, and, more generally, $X[k]=X[k+N]$.


Understanding multirate filtering fundamentals I am having trouble grasping some of the fundamental concepts of multirate filtering. I see from various sources that the basic building blocks of a multirate filter are the dyadic analysis and synthesis blocks.

Question 1:

The analysis block structure looks like the following, where the wide band signal is split into lowpass and highpass bands, each with a cutoff of FS/4 (Nyquist/2). Each band is then decimated by a factor of 2.

How can you accurately represent the signal in the high frequency band when it contains frequency information above the Nyquist limit of the new decimated sample rate?
Question 2:

The analysis block structure looks like the following, where the sub-band signal is interpolated, re-filtered and then summed.

What is the purpose of the second filtering?

I will answer question 2 first, and hopefully that will help explain what is going on with question 1.

When you sample a baseband signal there are implicit aliases of the baseband signal at all integer multiples of the sampling frequency, as shown in the picture below.

The solid image is the original baseband signal, and the aliases are represented by the dashed images.  I chose an assymetric (i.e. complex) signal to help demonstrate the inversion that happens at odd multiples of the sampling frequency.

You might ask, "Do the aliases really exist?"  It's a bit of a philosophical question.  Yes, in a mathematical sense they do exist, because all of the aliases (including the baseband signal) are indistinguishable from each other.

When you upsample by inserting zeros in between the original samples, you are effectively increasing the sampling rate by the upsample rate.  So if you upsample by a factor of two (putting one zero in between each sample), you increase your sampling rate and the Nyquist rate by a factor of 2, resulting in the picture below.

As you can see, one of the implicit aliases in the earlier image has now become explicit.  If you FFT the samples it will show up.  A non-rigorous proof that the DFT transform does not fundamentally change is given below.

Now that you have the two explicit aliases, if you just want the baseband alias then you have to low-pass filter to get rid of the other alias.  Sometimes, though, people use the other aliases to do their modulating for them.  In that case you would high-pass filter to get rid of the baseband signal.  I hope that answers question 2.

Question 1 is basically the inverse of question 2.  Suppose that you are already in the situation shown in the second picture.  There are two ways to get the baseband signal that you want.  The first way is to low-pass filter (thereby getting rid of the higher alias) and then decimating by a factor of two.  That gets you to picture #1.

The second way is to high-pass filter (getting rid of the baseband alias) and then decimating by a factor of two.  The reason that this works is that you are intentionally aliasing the signal into the baseband, thus, once again, getting you to picture #1.

Why would you want to do it that way?  Because in most situations the signals will not be the same, so you can choose which signal you want, or do them both separately.

If you are studying multi-rate processing, I highly recommend getting "Multirate Signal Processing for Communication Systems" by Frederic Harris.  He does a really good job of explaining the theory without neglecting the math, and giving a lot of practical advice too.

EDIT: Intentionally sampling a signal at less than the Nyquist rate is called undersampling.  The following is my attempt at explaining mathematically why the FFT does not change when you upsample.  "x[n]" is the original set of samples, "u" is the upsampling factor, and "x'[n]" is the upsampled set of samples.

$$
\begin{eqnarray*}
X[k] &=& \sum_{n=0}^{N-1} x[n]e^{-i2\pi kn/N}\\
X'[k] &=& \sum_{n=0}^{uN-1} x'[n]e^{-i2\pi kn/uN} , \{ &x&'[n] = 0, n \neq mu, m \in (0..N-1)\\
&x&'[n] = x[n/u], n = mu\\
&=& \sum_{n=0}^{N-1} x'[un]e^{-i2\pi kun/uN}\\
&=& \sum_{n=0}^{N-1} x[n]e^{-i2\pi kn/N}\\
&=&X[k]
\end{eqnarray*}
$$

Apologies for the ugly formatting.  I am a LaTex noob.

EDIT 2: I should have pointed out that the DFT's of x[n] and x'[n] are not truly identical.  The sample rate is higher, which as I explained in the earlier part of the answer, causes aliases to be "exposed".  I was trying to point out in my non-mathematician way that the DFT's are, aside from the sample rate, the same.


Concept of combining multiple FIR Filters into 1 FIR filter I try to combine my 7 band FIR filters into 1 FIR in order to process an audio signal.  I  use convolution to combine their coefficients but the effective is not the same as I apply 7 different FIR filter.  I also have the gain problem

H = conv(conv(conv(conv(conv(conv(h1,h2),h3),h4),h5),h6),h7)

My questions are
Is convolution a correct way to do combine multiple FIR filters into one filter
If not, then what technique should I look into?
If yes, then what did I do wrong?

Regards

What you tried to do in your example is combine the filters in series rather than parallel.
That is, if your signal is $$y = ((x*h_1)*h_2)\ldots*h_7$$ then, according to the associativity property of the Convolution, $$y=x*\underset{H}{\underbrace{(h_1*h_2\ldots*h_7)}}$$

However, if you wish to perform a parallel filtering then your signal is $$y=\sum_{i=1}^7 x*h_i$$
and acording to the distributivity property of the Convolution, $$y=x*\underset{H}{\underbrace{\sum_{i=1}^7 h_i}}$$


How do we know that CTFT of the autocorrelation function is the PSD? I know that the Fourier Transform of the autocorrelation function is
the Power Spectral Density. But how can we arrive at such a result intuitively?
Is it just a theorem?

$$\begin{align}
R_x(\tau) &=\int_{-\infty}^\infty x(t)x^*(t-\tau)\,\mathrm dt\\
&= \int_{-\infty}^\infty \left[\int_{-\infty}^\infty X(f)e^{j2\pi ft}\,\mathrm df\right]x^*(t-\tau)\,\mathrm dt\\
&= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x^*(t-\tau) e^{j2\pi ft}\,\mathrm dt\right]\,\mathrm df\\
&= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x^*(t-\tau) e^{j2\pi f(t-\tau)}\,\mathrm dt\right]e^{j2\pi f\tau}\,\mathrm df\\
&= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x^*(\lambda) e^{j2\pi f\lambda}\,\mathrm d\lambda\right]e^{j2\pi f\tau}\,\mathrm df\\
&= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x(\lambda) e^{-j2\pi f\lambda}\,\mathrm d\lambda\right]^*e^{j2\pi f\tau}\,\mathrm df\\
&= \int_{-\infty}^\infty X(f)X^*(f) e^{j2\pi f\tau}\,\mathrm df\\
&= \int_{-\infty}^\infty |X(f)|^2 e^{j2\pi f\tau}\,\mathrm df
\end{align}$$


Frequency shifting of a quadrature mixed signal I've a quadrature mixed signal with 10 MHz as highest frequency available. I want to shift it to any frequency lesser than 10 MHz. I should get it by multiplying this signal by cosine of the required frequency. Right?

If you have a quadrature mixed signal, you have two outputs of the signal in 90 degree phase to each other. To complete a full complex frequency conversion (which would shift the frequency in one direction, as opposed to a real frequency conversion which would create two sidebands) you need to multiply one of the outputs by the cosine term of the required frequency shift and the other output by the sine of the required frequency shift. By swapping sine and cosine ports, you change the direction of the shift.  The outputs after each multiplier are then summed in phase. See diagram and derivation using trig identities for the multiplication of sines and cosines below.

If you are comfortable with the idea of exponential frequencies (full complex frequency including positive and negative frequency etc) then I prefer the view point below to describe single sideband frequency translation. Comparing the two diagrams shows the value of converting "I" and "Q" signals as typically constructed in hardware to implement complex signals into one complex signal path for the purpose of analysis. Using cosines and sines instead of exp can be quite cumbersome especially as the complexity of the system increases.

For those less familiar with the notation, $Ke^{j\theta}$ is simply a complex value with magnitude K and angle $\theta$. Therefore $e^{j\omega t}$ is a signal whose magnitude is 1 and angle is increasing linearly with time at rate $\omega$; a spinning phasor counterclockwise on a complex IQ (polar) plot, representing a "positive" frequency. Similarly $e^{-j\omega t}$ spins clockwise and represents a "negative" frequency. With that, one should be well equipped to use the more compact and analytically simpler exponential notation.

Notice too how the two block diagrams that I included are related, with the second one taking the real of the product to make the outputs equivalent. To do frequency translation, we use complex conjugate multiplication (the conjugate assures we shift in the direction desired). Consider a full complex conjugate multiplication, as is done prior to the real operation:

$ e^{j\omega_ct} = I_1+jQ_1$     representing the carrier after the quadrature split

$ e^{-j\omega_st} = I_2-jQ_2$     representing the conjugate of the I and Q terms of your signal

$ e^{j\omega_ct}e^{-j\omega_st} = (I_1+jQ_1)(I_2-jQ_2) = (I_1I_2+Q_1Q_2)+j(I_2Q_1-I_1Q_2)$

The real portion is also called the dot product and the imaginary portion is called the cross product. Notice the dot product does match the implementation block diagram shown at the top of this answer.

Further details for the very interested to understand complex frequency translation:

The following additional figures and content were added to help answer the question on the subsequent down-conversion and to build a better understanding of the complex (positive and negative) frequency domain in describing frequency translation implementations.

First consider a generic quadrature mixed signal, it is at "baseband", meaning not modulated to a carrier for purpose of transmission, and centered about DC (0 Hz) with the positive and negative frequencies independent of each other. The way we implement such a signal in practice is with two real signals (we could choose the real signals to have one represent the real and the other the imaginary such as I and Q in I+jQ, or one represent the magnitude K and the other represent the phase $\theta$ such as $Ke^{j\theta}$. In any event, by observing the spectrum and noting that the positive and negative frequencies do not match, we know from that detail it is a complex signal and two real signal paths are required to represent this in implementation. For a real signal in contrast, the positive and negative frequencies would be conjugate symmetric, meaning identical in magnitude and opposite in phase.

Frequency spectrum of baseband quadrature modulated signal (this could be a single tone on one side, or multiple independent frequencies, point is the positive and negative frequencies are independent):

Now consider a carrier signal $cos(\omega_ct)$, this has an impulse in the frequency domain at the positive frequency $\omega_c$, and another at the negative frequency $-\omega_c$ (also shown as two complex frequencies using Euler's identity $2cos(\omega_ct) = e^{j.\omega_ct}+e^{-j\omega_ct}$):

If we split the carrier in quadrature (with a Hilbert transformer, or a 90° splitter, or just generate a sine and cosine component, but in the end have two tones with a 90° relationship to each other), then in implementation we have two real signals used to represent a single complex frequency $e^{j\omega_c}$  (you may start to see how it is so much easier doing all analysis with the e's and then do the implementation with the sines and cosines):

UPCONVERSION

To up-convert the baseband signal to the carrier frequency, we multiply the two signals in the time domain (which is convolution in the frequency domain; so for an impulse function such as our carrier it is a simple shift). First consider what would happen if we multiplied the baseband signal with the carrier $cos(\omega_ct)$ directly (eliminate the 90° splitter in my first figure, replacing $sin(\omega_ct)$ with $cos(\omega_ct)$). Notice the output would still be complex as the positive and negative frequencies are not the same- so we could not transmit this with a single antenna (2 real signals are required to represent a complex signal). If we proceeded to take the real portion,
distortion of our signal would result as depicted in the figure, so we can't transmit that either.

Figure showing spectrum result with an INCORRECT upconversion approach:

This explains why we need to also take the Hilbert transform of the carrier (converting the cosine into sine and cosine components, and representing a  complex version of the carrier that is a single spinning phasor $e^{j\omega_ct}$, or in the frequency domain a single impulse in the positive domain with no negative component) prior to multiplication. We take the real portion of the product and can then send an undistorted real waveform representing our baseband signal modulated up to the carrier frequency.

Figure showing spectrum result with a CORRECT upconversion approach:

DOWN-CONVERSION

With the details of the up-conversion process in mind, it will now be very easy to answer the question if a quadrature conversion of the received signal is required prior to multiplying by sine and cosine to recover the baseband signal. The answer is no, but there is an advantage to doing so in situations when an image frequency cannot be easily filtered.

Figure showing spectrum doing a downconversion with no Hilbert transform of the received signal. It is interesting to note if you swapped sine and cosine in the (real) multipliers, this would represent a single impulse in the positive frequency domain, shifting the negative frequency component of the received signal to baseband, causing the spectrum at baseband to be inverted:

Figure showing spectrum doing a downconversion with a Hilbert transform of the received signal prior to sine/cosine multiplication. Here you can see what would happen if you swapped an input to one of the (real) multipliers: The signal would be shifted to a higher frequency twice the carrier instead of to baseband.


DTMF Tone Generation in MATLAB not sounding right I'm starting out with some DTMF Tone Generation using Matlab and I wanted to find some theory behind how to generate the various tones for a 4x3 keypad consisting of 0-9, * and #. So far I have the following function:

function y = OutputTone(c)

t = 0:0.1:500;
num = double(c);
%fs = 8000;

if (num == 49) || (num == 50) || (num == 51) % 1, 2, 3
    f1 = 697;
    if (num == 49)
        f2 = 1209;
    elseif (num == 50)
        f2 = 1336;
    elseif (num == 51)
        f2 = 1477;
    end
elseif (num == 52) || (num == 53) || (num == 54) % 4, 5, 6
    f1 = 770;
    if (num == 52)
        f2 = 1209;
    elseif (num == 53)
        f2 = 1336;
    elseif (num == 54)
        f2 = 1477;
    end
elseif (num == 55) || (num == 56) || (num == 57) % 7, 8, 9
    f1 = 852;
    if (num == 55)
        f2 = 1209;
    elseif (num == 56)
        f2 = 1336;
    elseif (num == 57)
        f2 = 1477;
    end
elseif (num == 42) || (num == 48) || (num == 35) % *, 0, #
    f1 = 941;
    if (num == 42)
        f2 = 1209;
    elseif (num == 48)
        f2 = 1336;
    elseif (num == 35)
        f2 = 1477;
    end
end

y = sin(2*pi*f1*t) + sin(2*pi*f2*t); % DTMF Signal in Time Domain

end

As the code describes above, I take in an input character (0-9,*,#), convert it to ASCII and then select the appropriate DTMF frequency pertaining to the key that is pressed. For example, if 1 is input to the function, then the frequencies f1 = 697Hz and f2 = 1209Hz will be selected and applied to the output sinusoid y.

However, when I play back the sound of y, it doesn't sound anything like the key tones I hear on my smartphone or in this link. What is the correct way to represent y?

You have commented the line that you assign fs. Then the time duration t should be defined with respect to fs.

It should be as follows:

fs = 8000;
t = 0:1/fs:500;


Relation between Bandwidth and Baud-Rate for 8-PSK I have a query on relationship between Baud Rate and Bandwidth for 8-PSK.

In 8-PSK as 3 bits together form a symbol, so Bandwidth = 2R/3 , R is Bit Rate

As bit rate is 3 times BAUD RATE for 8-PSK, hence Bandwidth has to be = 2* (BAUD RATE)

But in many places is see the below statement ::::

For PSK the baud rate is the same as the bandwidth I was going through this - This

How can this be possible.??

The final statement is correct, for PSK with proper pulse shaping the baud rate and the bandwidth are the same (the bandwidth will typically be 20-30% higher than the symbol rate, but read on). The baud rate is the symbol rate or the rate of change from one symbol to the next in the modulated constellation.

Your earlier figures suggest unfiltered waveforms and "bandwidth" definition to be the first main lobe of the resulting Sinc function spectral shape. Specifically, with a random data pattern, the spectrum for an unfiltered PSK waveform will be a sinc function with a null to null bandwidth from the first null that is twice the baud rate.  (and power spectral density will be sinc squared). Perhaps this is the source of your confusion.

Once filtered (typically with raised cosine or root-raised cosine filters or other similar Nyquist filters) the bandwidth can approach the symbol rate, but will be slightly wider given the parameters of the Nyquist filter. (For an infinitely long filter the bandwidth can be exactly the symbol rate but that is not practical).

This also applies for all the same reasons for QAM which is shown in the figure below (I happen to have figures for QAM that help demonstrate this), and when you see and understand the additional explanation below you will know how it can be applied to determine the bandwidth of other modulated waveform. The figure below shows a QAM waveform versus time (the real portion of the waveform) as an unfiltered waveform in blue versus what it would look like with raised cosine pulse shaping in red, and the resulting spectrum for each case is shown directly below that. This waveform was modulated with a baud rate of 1 MS/s.

The unfiltered view is clear from observing the resultant modulated signal as a stream of impulses that are the convolved with a rectangular pulse that is the width of your symbol period. Again see the example related to QAM below, which equally applies to PSK. If the waveform was not pulse-shaped, the stream of impulses representing the modulation symbols would be convolved in time with a rectangular function of width $T$ (where $T$ is the symbol period) to result in the blue waveform above. Note for PSK this stream of impulses is complex values all with same magnitude by varying phase from symbol to symbol. The stream of impulses if indeed completely random and independent from symbol to symbol would have a white spectral density. Since the impulses are convolved with a rectangular function in time, its spectrum is multiplied by the Fourier transform in frequency, which for a rectangular function is a Sinc with a first null at $1/T$! Thus it should be very clear that the spectrum for a unfiltered PSK or QAM modulated waveform will be a Sinc function with a main lobe of twice the symbol rate. Pulse shaping reduces the bandwidth further.

For an immediate insight into why the bandwidth can reduce to half the main lobe at most, consider a modulation case where the waveform alternates on every symbol between two symbols. The result would be a square wave at a rate of half the symbol rate. This would represent the highest frequency data pattern possible (such as for BPSK the data pattern of 101010 vs 11001100), and filtering down to half the main lobe would allow the fundamental of this square wave to still pass, as a sine wave yet still containing the information (101010) for subsequent demodulation. Filtering any lower than half the main lobe and this particular pattern (representing the highest frequency possible) would be completely removed.


FFTW - Inverse FFT input data I am trying to use FFTW function for IFFT:

fftw_plan_dft_c2r_1d(size, out_cpx, out, FFTW_ESTIMATE).

My goal is to input desired frequency response data and obtain related coefficients using Inverse Fourier Transform. Since out_cpx is complex, I calculated magnitude (0,1) and phase, taking into account it is linear phase filter. Since I know the phase I calculated real and imaginary parts as Acos(w) and Asin(w) and with these pairs I populated out_cpx[i][0] and out_cpx[i][1]. When I perform Inverse Fourier Transform using fftw_plan_dft_c2r_1d() I assume I am right away getting FIR coefficients h[n] as the result. With these coefficients i go to MATLAB and get response using freqz(C,1) where C holds my coefficients. But I am getting a lot of ripples in my graph for response and for phase. Am I doing something wrong here? It is not clear to me how to specify phase of the desired frequency response.

The ripple you see is due to Gibbs phenomenon and is caused by your filter having a step response. You will have to pass the frequency response through a window function first.


Using IFFT for the purpose of "Multiple Orthogonal Frequency Shift Keying" Is it possible to create a set of complex orthogonal sinusoids using FFT or IFFT
In other words a kind of  "Multiple Frequency Shift Keying" using  multiple sinusoids that are created with a fourier transform ( FFT or IFFT) and where one frequency is sent at a time?

Yes, basically this is OFDM with an on/off constellation for each subcarrier instead of a QAM constellation as often used. So you would map your data to FFT values (each FFT bin is either 1 or 0), then take the IFFT to create the transmit waveform, with the cyclic prefix etc as in standard OFDM implementations.


How to do an FFT of a line profile obtained by imagej? First of all, I'm a total "signal processing" newbie. So, bear with me.

I know that in ImageJ/Fiji I can use the Plot Profile function to select a line profile of an image. Is there a simple way to plot a 1D FFT transform of the data along the line?

If there is a better SO site for this kind of question, please let me know.

It appears there is no straightforward way to do it from the UI. After looking through ImageJ docs and finding out how to start the interpreter, I've come up with something that is essentially this:

macro "FFT1D Tool - C099A8877099L0088" {
    profile = getProfile();
    fft=Array.fourier(profile,"None");
    s=newArray(lengthOf(fft));
    for (i=1;i<lengthOf(fft);i++) {
        s[i]=2*PI/i;
    }
    s[0]=0
    Plot.create("Fourier spectrum of the line profile","wavenumber","FFT(line)",s,fft);
    Plot.show();
}

If placed in StartupMacros.fiji.ijm in the fiji folder, on the next start of fiji, the FFT1D tool will show up in a toolbar menu and, after a line was drawn over the image, a click on the toolbar button will execute the macro.

A caveat: DFT works great for periodic functions and input line sizes that are powers of two. Manually drawing over the image is likely going to result in a non-power-of-two, aperiodic array. Running this through a DFT will result in Gibbs effect which will ruin the accuracy of the transform.    Window functions (substitute "None" for one of the available window functions) can help ameliorate this effect. It helps if the intensity of the image decays towards the endpoints of the profile line.

Also, check this example of using the Array.fourier function.


Correlation of two signals at different sample rate I have a signal $X$ and I want to find signal $Y$ in $X$, I can achieve this by cross-correlation. Signal $X$ is sampled at $8000\textrm{ Hz}$, and $Y$ signal is sampled at $44100\textrm{ Hz}$.

My question is "should both signals be sampled at same rate to correlate to find $Y$ in $X$ " Or is it okay to correlate two signals sampled at different rate to find  $Y$ in $X$?

If you are looking for the exact signal, then yes you need to resample one of them to the same rate as the other before performing a cross correlation. Otherwise you will be cross correlating against a time-stretched (or compressed) version of the desired signal (which probably will not yield accurate results). If your $Y$ signal bandwidth is less than $4000\textrm{ Hz}$, then you may want to use an anti-aliasing filter and resample it to $8000\textrm{ Hz}$ for reduced computational load in the cross correlation. If not, then you must resample the signal $X$ to $44100\textrm{ Hz}$ in order for a cross correlation to be meaningful.


Implementation of Block Orthogonal Matching Pursuit (BOMP) Algorithm - Fix Given Code This is my implementation which doesn't work:

function [result,gamma,T]=block_OMP4(phi,signal,N,d,maxIter)
[m n]=size(phi);
r=signal;
gamma=[];
T=[];
i=1;
tol=1e-4;
err=1;
x=zeros(N*d,1);
val=[];
figure,
while  (i<=maxIter && err>tol)

    maxv=[];
    for j=1:N
        startx=(j-1)*d+1;
        endx=startx+d-1;
        block=phi(:,startx:endx);
        t=block'*r;
        maxv(j)=norm(t,2);
%
    end
    [maxvl,maxvindx]=max(maxv);
    gamma=union(gamma,maxvindx);
    x2=lsqr(phi,r);
    xb=zeros(N*d,1);
    for i=1:length(gamma)
        stratx=(gamma(i)-1)*d;
        endx=stratx+d-1;
        xb(stratx:endx)=x2(stratx:endx);
        x(stratx:endx)=x(stratx:endx)+xb(stratx:endx);
    end

    y=phi*x;
    r=signal-y;

    val(i)=norm(r);
    err=norm(r)
    i=i+1;

end
result=x;
end

Anyone has a working implementation?

Answer taken from Implementation of Block Orthogonal Matching Pursuit (BOMP) Algorithm.

The Block Orthogonal Matching Pursuit (BOMP) Algorithm is basically the Orthogonal Matching Pursuit (OMP) Algorithm with single major difference - Instead of selecting single index which maximizes the correlation we chose a set of indices, sub set of columns of the matrix and the solution vector.

A good reference for the algorithm is given in:

An Optimal Condition for the Block Orthogonal Matching Pursuit Algorithm.
Block Sparsity: Coherence and Efficient Recovery.

The code is given by:

function [ vX ] = SolveLsL0Bomp( mA, vB, numBlocks, paramK, tolVal )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveLsL0Omp( mA, vB, paramK, tolVal )
% Minimizes Least Squares of Linear System with L0 Constraint Using
% Block Orthogonal Matching Pursuit (OMP) Method.
% \arg \min_{x} {\left\| A x - b \right\|}_{2}^{2} subject to {\left\| x
% \right\|}_{2, 0} \leq K
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix (Fat Matrix). Assumed to be
%                           normlaized. Namely norm(mA(:, ii)) = 1 for any
%                           ii.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numBlocks         -   Number of Blocks.
%                           The number of blocks in the problem structure.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - paramK            -   Parameter K.
%                           The L0 constraint parameter. Basically the
%                           maximal number of active blocks in the
%                           solution.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - tolVal            -   Tolerance Value.
%                           Tolerance value for equality of the Linear
%                           System.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range [0, inf).
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  An Optimal Condition for the Block Orthogonal Matching Pursuit
%       Algorithm - https://ieeexplore.ieee.org/document/8404118.
%   2.  Block Sparsity: Coherence and Efficient Recovery - https://ieeexplore.ieee.org/document/4960226.
% Remarks:
%   1.  The algorithm assumes 'mA' is normalized (Each column).
%   2.  The number of columns in matrix 'mA' must be an integer
%       multiplication of the number of blocks.
%   3.  For 'numBlocks = numColumns' (Equivalent of 'numElmBlock = 1') the
%       algorithm becomes the classic OMP.
% Known Issues:
%   1.  A
% TODO:
%   1.  Pre Process 'mA' by normalizing its columns.
% Release Notes:
%   -   1.0.000     19/08/2019
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

numRows = size(mA, 1);
numCols = size(mA, 2);

numElmBlock = numCols / numBlocks;
if(round(numElmBlock) ~= numElmBlock)
    error('Number of Blocks Doesn''t Match Size of Arrays');
end

vActiveIdx      = false([numCols, 1]);
vR              = vB;
vX              = zeros([numCols, 1]);
activeBlckIdx   = [];

for ii = 1:paramK

    maxCorr         = 0;

    for jj = 1:numBlocks
        vBlockIdx = (((jj - 1) * numElmBlock) + 1):(jj * numElmBlock);

        currCorr = abs(mA(:, vBlockIdx).' * vR);
        if(currCorr > maxCorr)
            activeBlckIdx = jj;
            maxCorr = currCorr;
        end
    end

    vBlockIdx = (((activeBlckIdx - 1) * numElmBlock) + 1):(activeBlckIdx * numElmBlock);
    vActiveIdx(vBlockIdx) = true();

    vX(vActiveIdx) = mA(:, vActiveIdx) \ vB;
    vR = vB - (mA(:, vActiveIdx) * vX(vActiveIdx));

    resNorm = norm(vR);

    if(resNorm < tolVal)
        break;
    end

end

end


Is there an equation that will graph the shape or predict the maximum amplitude of a single sample impulse fed through a One-Pole LPF? If you feed a single sample impulse of a given amplitude through a One-Pole LPF you get an envelope that looks like this:

The output of the LPF from its peak will settle into an exponential decay where time to reach $1/e$ amplitude is given by $1/(2πf)$ where $f$ is the cutoff frequency of the LPF.
I am wondering if there's any way to predict the maximum amplitude of the envelope from the sample rate, LPF cutoff freq, and amplitude of the single sample impulse. In general, is there a non-recursive equation that can describe or approximate the curve that results from these three factors?
If such an equation can be developed or exists, then one can take the derivative and solve where it equals zero to find the max.
Is this possible?
Alternatively, I can solve it numerically by just iterating through a temporary setup of a single sample impulse into an LPF of the given parameters until (output_1 < output). This will find the max as well but this is tedious and computationally consuming.
I would hope some equation can do this faster.
Thanks for any thoughts.

As mentioned in the comments, an impulse response with the shape shown in your question can only be obtained by a system with two real-valued poles. In continuous time, with two distinct poles, the total impulse response is
$$h(t)=\frac{e^{-\alpha_1t}-e^{-\alpha_2t}}{\alpha_2-\alpha_1}u(t)\tag{1}$$
It's a basic exercise to determine the location of the maximum of $(1)$, and the result is
$$t_{max}=\frac{\log\left(\frac{\alpha_1}{\alpha_2}\right)}{\alpha_1-\alpha_2},\qquad\alpha_1\neq\alpha_2\tag{2}$$
For $\alpha_1=\alpha_2=\alpha$, the total impulse response is
$$h(t)=te^{-\alpha t}u(t)\tag{3}$$
and the location of the maximum is
$$t_{max}=\frac{1}{\alpha}\tag{4}$$
The figure below shows an example of an impulse response $(1)$ with $\alpha_1=0.4$ and $\alpha_2=0.6$, and with the location of the maximum predicted by $(2)$.


How to plot noise shaped spectrum of First order Incremental Sigma Delta ADC's output? A little background:

This question is about viewing the output spectrum of incremental delta sigma modulators (IDSM). IDSM is a subcategory of general delta sigma modulators where the integrators are reset for every conversion cycle with the advantage of being multiplexable along several channels. (https://pdfs.semanticscholar.org/1f62/af5ec0ba627e8ef5ebf81a479307b059a116.pdf).

I am aware that a general sigma-delta noise shaping architecture has a noise shaping characteristic of a differentiator as shown below:

1,2 and 3 indicate the orders of noise shaping and sigma delta architectures too.

For example, if we consider the following continuous-time sigma-delta architecture (first order), the noise shaping is governed by the equation:

NTF(z)=(1-z^(-1)), which translates to the curve (1) of the first figure.

Simulation background for IDSM:

When I simulate an incremental sigma delta converter (IDSM), I run it for (N*(1/fs)*OSR), where it runs for OSR clock cyles for single conversion. I run the output through the counter and then scale it by (1/M) for every reset (Decimation) and then view the spectrum with N points, I am able to get the desired output SNR for an IDSM.

One conversion cycle is shown below (M=OSR).

Core question:

I designed a continuous-time first-order sigma delta modulator in matlab and the output spectrum is shown below:

Here I am able to see the first order noise shaping (-20dB/decade noise increment).

Now, how do I simulate the output from an IDSM to view the result of noise shaping?.

I saw some articles, plotting noise shaped spectrum for IDSM (in free running DSM mode). I am really not sure how to run an IDSM in this mode, because logically it would require a reset for each conversion and therefore I am stuck.

(Again, I could decimate the output of IDSM and then view the spectrum - which gives me the correct SNR. But I'm interested in viewing the noise shaping behavior.)

If what i am expecting is not achievable, please let me know as well - that would be helpful too.

Update: as Sundar commented his question is specific to the “incremental Sigma Delta ADC” of which I only have a vague understanding of. To the extent the incremental implementation is a traditional Sigma Delta ADC that is reset periodically I would believe this model still applies as described similar to using a moving average frequency response to model an integrate and dump. Comments welcome.

Below is how I use an equivalent loop equation to model the noise shaping in first and higher order Sigma Delta architectures. This is qualitative since the quantizer is highly non-linear (So would a apply from a small signal perspective, and can be further more accurately approximated with a piece-wise linear approximation to the quantizer gain.

Below shows continuous time equivalents, with $N_q$ as the quantization noise input from the 1 bit converter. The discrete time equivalent is easily derived by mapping the integrator. An integrator using either Backward Euler or Impulse Invariance mappings both result in:

$$\frac{1}{s} \leftrightarrow \frac{T}{1-z^{-1}}$$


Confusion understanding digital signal processor? I am reading wikipedia article about digital signal processor and as shown in end of attached snapshot
." Most general-purpose microprocessors can also execute digital signal processing algorithms successfully, but may not be able to keep up with such processing continuously in real-time"
Please elaborate with example/application? What is meant by real-time here in context of dsp?

The staple processor instruction of specialized DSP processors is a multiply-accumulate instruction. The vector multiplication and summation are performed in a single command fetch, it's a sort of SIMD (but not exactly). Only advanced MCU's have this instruction in their command set. Hence the advantage of DSP processors over general-purpose MCU's w.r.t. ALU performance on typical signal-processing operations, even when comparable technological process are used in their manufacturing. The same for FPGA devices with internal DSP blocks. In short, less command fetch cycles, higher performance.


Fuse two sources of linear acceleration with a Kalman filter How would I fuse two different sources of linear acceleration with a Kalman filter (perhaps linear acceleration readings from an IMU and from a dedicated accelerometer)?

My state is defined by distance and velocity; it is very similar to the following post, except for my availability of two sources of linear acceleration:

Using the Kalman filter given acceleration to estimate position and velocity

I'm not sure which acceleration reading I would use for my control vector ($u_{k}$) in the predict function.

That depends on how many of the accelerometer parameters (mostly drift and misalignment) you're trying to estimate.

If the IMU and the 'extra' accelerometer were in perfect alignment (and if their statistics are Gaussian), then the optimal combination of their outputs would be a simple weighted sum: $\vec {\hat a} = k_1 \vec a_1 + k_2 \vec a_2$ where (assuming equal scaling), $k_1 = \frac{\sigma_2^2}{\sigma_1^2 + \sigma_2^2}$, $k_2 = 1 - k_1$, and $\sigma_1$ and $\sigma_2$ are the noise variance of each accelerometer.

If you look at that first equation, you'll see that if one accelerometer is markedly quieter than the other, then unless it has other problems (misalignment, scaling, initial offset, etc.) your best bet is to just ignore the noisy accelerometer altogether.  Even at the point where one accelerometer has three times less variance than the other, the ultimate contribution of the "noisy" accelerometer ends up being $1/9^{th}$ of the total -- which is definitely tending into "why bother?" territory.

If you've got some oddball combination, i.e. you've got a really quiet accelerometer with terrible initial offset or misalignment vs. a really noisy one with low offset and misalignment, then life gets more complicated -- in fact, I'd have to ponder on exactly how you'd express that in a way that would make a sensible Kalman filter.


Accuracy of IFFT(FFT(X)) unit transformation I want to use FFTs to do seriously accurate interpolation on band-limited data.  To do that, I need to get a handle on the fundamental accuracy of the FFT() and IFFT() algorithms available.  My idea is to look at IFFT(FFT(X)) which should be a unit transformation.  I want to get the maximum difference between the input and output data down below the region of 1 part in 10^8 or even 10^9.  Currently I'm achieving ~1 part in 10^6 using 64-bit float, but I don't know if that is in line with what is to be expected.  I'm thinking that surely someone has done all this before.  Does anybody know the answers?  What can expect if I try using higher-precision data types?  With 64-bit float speed is not a concern.

Your mileage may vary, but in multiple environments on my x86_64 bit linux machine, I get much less error than 1 / 10^9.

julia:

julia> a=randn(32768); f=ifft(fft(a)); print(maximum(abs(f-a)))
1.782802105510316e-15

octave:

octave:1> a=randn(32768,1); f=ifft(fft(a)); disp(max(abs(f-a)))
1.8184e-15

numpy:

In [1]: %pylab
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

In [2]: a=randn(32768); f=fft.ifft(fft.fft(a)); print(max(abs(f-a)))
7.58800755928e-15

The default for all these environments is float64.  If you want to share how the signal is bandlimited, that might shed more light on the problem.

UPDATE: To satisfy my curiosity, I tried this with a single sine-wav and basically get the same results:

julia> th=collect(0:32767)*(1000*pi/32767); a=sin(th);
       f=ifft(fft(a)); print(maximum(abs(f-a)))
9.393934579042684e-16

Again, without more knowledge about the signal or the computation environment, this is hard to address.


Incoherence: Compressed Sensing (CS) vs Matrix Completion (MC) I am seeking a clarification of the concept of Incoherence within the MC framework.

Specifically,

1) the literature mentions the application of a "strong incoherence" given a set of assumptions. Intuitively, how is the strong incoherence different than "incoherence".

2) Does the Restricted Isometry Property have anything to do with low rank Matrix Completion of natural images, or does it only pertain to the reconstruction of signals in CS (i.e. y = Ax).

Background

There is an incoherence property $\mu$ used in CS to measure the correlation between the sensing and representation bases; whereby low correlation enables signal reconstruction of sparse signals with fewer samples and high correlation makes it nearly impossible regardless of the sparsity of the signal.

For MC, there is also an incoherence property $\mu$: it measures the correlation between the singular vectors of the matrix (U and V from SVD) and the standard basis vectors (i.e. $\ e_i$ = vector of zeros except for a 1 at index i). As with CS, low rank (i.e. singular value-sparse) matrices exhibiting a low correlation are more likely to be reconstructed than those with high correlation.

Both of these incoherence properties are used to stipulate the minimum number of samples required to reconstruct a signal (1D, 2D, respectively). However, I cannot see where the matrix incoherence property addresses the nature of the sampling; I.e. it only tells me which matrices are more likely to be reconstructed. I believe this is because the "nature of the sampling" that is assumed (by CS community doing MC) for low rank MC problems is strictly "uniform random sampling".

Now, restricting our attention only to those low rank matrices with low $\mu$, is there an additional MC incoherence property that characterizes the nature of the sampling to say whether or not a matrix can be reconstructed? (I'm looking for a link back to CS where they compare the relationships between different representation/sensing bases to see if they work well together to reconstruct the signal.)

Thank you.

Partial answer:

2) Matrix completion does "not rely on any notions of restricted isometry." Candes and Recht, 2008.


Adaptive Predictive Coding: Transmit signal AND prediction coefficients I'm not well versed in the literature so please forgive/inform me if my terminology is wrong:

I have a bounded real-valued, discrete-time signal $x_n \in \mathbb [0, 1]$, and I want to (lossily) transmit this signal via a bitstream $q_n \in \{0, 1\}$, so that it can be reconstructed into $\hat x_n$, and the average reconstruction error $\frac1N \sum_{n=1}^N (x_n - \hat x_n)^2$ is minimised.

$x_n$ will in general be temporally redundant, and I'd like to exploit this redundancy to make the most efficient use of my bits.  In my case, I expect $x_n$ to settle towards a fixed point for a long time and then suddenly jump towards a new fixed point (as you may expect from the trace of a single pixel in video).

Example Signal: It goes though "settling" phases where it remains fairly static, and "dynamic" phases where it transitions quickly towards a new settling point.

From what I gather, Predictive Coding is the standard approach here - subtract the predictable component of the signal, encode the rest, and reconstruct the predictable component on the other end.  The trouble is that the optimal predictive coding parameters (at least for a Linear Predictive Coder) will change depending on the state of the signal: When $x$ is dynamic, I expect to encode it at a coarse resolution.  When $x$ is static, I expect to use my bits to communicate fine increments to the value of $x$.  So I need to communicate my predictive coding coefficients along the same bitstream $q_n$.

Is there a standard approach for this kind of problem?

This sounds like a variant of Differential Pulse-Code Modulation. In this set of techniques, the correlation of adjacent values and the effect of quantisation are exploited for "compression" (whether this is for transmission or storage). That is, instead of encoding the sample values, you encode the difference between them which, within limits, could be described by a much smaller number than the values themselves.

However, another technique that might be useful to you here is run length encoding. In run length encoding, the encoded word describes either the sample value itself or how many repetitions of the same value are to be decoded.

In either case, you would require higher bandwidth around the signal transients (many messages sent between the encoder and the decoder) and much lower during the longer times of constant values.

Hope this helps.


2D Convolution in the Spatial Domain vs Frequency Domain Suppose, I have this kernel.
-1, -1, -1,
-1,  9, -1,
-1, -1, -1

Can this kernel be used in a FFT based convolution? How?
What could be the reason of my failure?

Related:

Image Convolution in Frequency Domain.
FFT Convolution - 3x3 Kernel

To answer your question:

Can easily be done.
One must remember that the short signal (The Kernel) must be padded (With zeros) to have the same size as the image before the DFT conversion.
Once they have the same size all needed is to convert into the Frequency Domain have element by element multiplication and transform back.
A side note would be that this way you assume periodic boundary condition.
I'm not a C# / C Coder.
I can provide you a MATLAB reference if needed.
Just make sure you do the following steps:

Pad the kernel to the image size (Pad it with zeros to the right and left).
Convert into the Frequency Domain.
Multiply Element by element.
Convert back to the spatial domain.

2 Remarks:

In cases of large image relative to the size of the kernel you better (Efficiency wise) apply it in the spatial domain.
The method above describes to do Circular Convolution (See Applying Image Filtering (Circular Convolution) in Frequency Domain). In order to apply Linear Convolution See How Much Zero Padding Do We Need to Perform Filtering in the Fourier Domain?

Enjoy...


How to calculate resolution of DFT with Hamming/Hann window? The frequency resolution of a DFT with a rectangular window of size $N$ is given by $f_s/N$. However, when using other window functions like a Hamming or Hanning window the resolution gets worse. How can I calculate the frequency resolution of a DFT when using a Hamming/Hanning window of size $N$?

I found some numbers regarding "effective noise bandwidth (ENBW)" but am not sure what exactly it means. The numbers were 1 for a rectangular window, 1.5 for Hanning and 1.37 for Hamming window. Does that mean the frequency resolution using a Hanning window is $1.5 \cdot f_s/N$?

The frequency resolution of a DFT is indeed significantly affected by the choice of window. Two key parameters to consider when selecting a window are the frequency resolution and dynamic range. The frequency resolution is the ability to discern two closely spaced frequencies at similar power levels (therefore a metric of what that spacing is would be the frequency resolution), while dynamic range is the ability to discern two signals at different frequencies and significantly different power levels; specifically the ability to see a low level signal in the presence of a much stronger signal (therefore a metric of what this power difference is would be the dynamic range). Typically a higher dynamic range is achieved by using appropriate windowing at the expense of frequency resolution. The highest frequency resolution that can be achieved is with a rectangular window (and the lowest loss; any other window will also introduce a processing loss).

The Equivalent Noise Bandwidth (ENBW) is a commonly used metric for frequency resolution. The ENBW is the bandwidth (often given in FFT bins) of a brickwall filter that would result in the same noise power as the DFT "filter" for white noise (the DFT can be described as a bank of filters). For the rectangular window, the ENBW is one bin, but the poorest dynamic range (sidelobes with peaks rolling off at a $1/f$ rate in magnitude.  Any other windowing used would have a higher ENBW.

For a given length N window function $w[n]$, the ENBW in bins is computed as follows:

$$ENBW = N \frac{\sum \left(w[n]\right)^2 }{\left(\sum w[n] \right)^2}$$

See this classic paper by fred harris tabulating the ENBW as well as other paramters of many commonly used windows. fred harris windowing

Below is a slide I have further depicting the ENBW for a rectangular window:

Below compares a N = 30 rectangular window (denoted the Dirichlet kernel) to a Blackman window of the same length. The decrease in level shown of -7.83 dB is the coherent gain for this particular window used. The coherent gain is simply the average of the window weights. The result below is is consistent with the coherent gain of 0.42 as listed for the Blackman window tabulated in the fred harris reference given; $20Log_{10}(0.42) = -7.54$ dB.  What is very clear from this picture is the significant increase in both dynamic range and frequency resolution.

The reason for not being an exact match to the tabulated results is the window length. As the window length approaches infinity the physical result approaches the tabulated, as shown in the figure below for coherent gain. For this reason caution must be applied when using the tabulated results for any smaller length windows.


Mixed time/frequency domain data I'm afraid this might be a bit of an elementary question, but I'm afraid I'm a little new to this. I've got a system that measures a frequency spectrum to find resonant frequencies. It does this by applying an impulse, and then performs a linear sweep through a range of frequencies measuring the amplitude of the response.

As far as I can tell, this seems to produce data that contains both frequency- and time-domain information, because the frequency being measured changes over time. The results seem to agree with this: the initial impulse decays over "time", with individual excitations being visible at various frequencies. See here for an example (frequency on the horizontal axis, amplitude on the vertical):

I would like to remove from this data the "time" component, leaving just the resonant responses. I have tried mirroring the data, then deconvolving a "top hat" function, and FTing the result, but I am not convinced this corresponds to anything meaningful. Could anyone provide any insight?

this is John BG

1.- You have to work with Signal Power, not just input or output:

Circuit resonance is measured with the quality factor Q and does not only have to do with the signal you send in, or the signal that bounces back, but with those AND the impedance that such signals find.

Literature reference Microwave Engineering (Pozar) 1 defines Circuit Resonance the following way

2.- Measure the impedance:

Without the impedance, there's no way you can assess Q.

You are focusing your effort on the test signal and it's spectrum, but there's not a single reference to the impedance of the circuit you are measuring.

There's no way you can measure Q without at least an additional signal, than can be the impedance that the test input faces, or alternatively the reflected signal, that again you haven't mentioned.

3.- Narrow band or wide band? doesn't really matter

Despite the definition of resonance is done with a central frequency f0, it may be the case that the circuit resonates on each of the frequencies within the Chirp signal you use to poke the system.

4.- Use Chirp with constant amplitude

From the amplitude of the signal you shown in the question, the amplitude decreases with increasing frequency, that complicates calculations.

5.- Wrapping Suggestion:

Why don't you

start with a constant amplitude input chirp,
measure the return signal
from there get the impedance
calculate the 'sharpness' of the impedance notch over frequency,
which is actually what defines how resonant is a circuit,
from there you can tell whether the circuit under test is resonant
throughout the chirp signal frequency band or only on certain
frequencies.
Once you have all this working, you may want to start playing with
the amplitude of the chirp input signal.

If you follow the above steps, please let us know the results you obtain updating this QA with further comments, whether it worked for you, or if not what code have you written on such regard.

If you send me by email the MATLAB code you used to generate the question signal, I will have a look and look for a way to get the reflected signal.

John BG

Literature Reference (1):

https://www.amazon.co.uk/Microwave-Engineering-Written-Publisher-Hardcover/dp/B00SLRKBV2/ref=sr_1_1?ie=UTF8&qid=1529281040&sr=8-1&keywords=David+Pozar+Microwave+Engineering

Excerpts of the book can be read here: http://www2.electron.frba.utn.edu.ar/~jcecconi/Bibliografia/Ocultos/Libros/Microwave_Engineering_David_M_Pozar_4ed_Wiley_2012.pdf

There's also a solutions manual available here: https://www.scribd.com/doc/176505749/Microwave-engineering-pozar-4th-Ed-solutions-manual


Cyclic convolution between two matrices in matlab How can I perform cyclic (periodic) convolution between two matrices in Matlab? And also, what would be the normalization term?

Assuming that cconv2 does not exist, you can do it with the Fourier transform and multiplication:

A = some_MxN_matrix;
B = some_KxP_matrix;
Mfft = max(M,K);
Nfft = max(N,P);
C = ifft2( fft2(A,Mfft,Nfft) .* fft2(B,Mfft,Nfft) );

Hopefully the above helps.


Downsampling and Convolution: When Can one Swap Between Them? (Code snippet with example below).

In MATLAB, assume I have a long signal vector y (length N) that I wish to convolve with a square wave h (consisting of H 1's). This would give a convolution result of length N+H-1.

Since I have to perform this operation quite often, I thought that I could speed up my code by downsampling the signal y (and also the square wave h) and then perform the convolution. The downsampled signals are referred to as y_ds and h_ds.

However, the result is shifted when compared to the downsampled version of the convolution result (if corrected for scaling differences). So in general I found that:

  conv(y_ds,h_ds) $\neq$ downsample(conv(y,h))

However, the shift is between the two results is not constant, i.e. shifting such that the maxima align, causes the edges to be different.

My questions

Where does this shift originate from?
How can this code be modified such that the downsampled convolution result equals the convolution result of the downsampled signals?

Kind regards

N = 100;
dsfactor = 3;
H = dsfactor*3;
x = linspace(0,10,N);
y = sind(18*x);
h = ones(H,1);

convyh = conv(y,h);

convyh_ds = downsample(convyh,dsfactor)./H;
convyh_ds3 = conv(downsample(y,dsfactor),downsample(h,dsfactor))./(ceil(H/dsfactor));

figure;
plot(convyh_ds,'DisplayName','downsample(conv)','LineWidth',1);
hold on;
plot(convyh_ds3,'DisplayName','conv(downsample)','LineWidth',1);
grid minor
legend('show');
line([0 length(convyh_ds2)],[0 0],'LineStyle','--');

I thought this was because you were getting aliasing in your signal during one of the downsampling operations. However, I re-wrote your code in R (included below) and I did not see the same effect  I suspect conv is doing something odd. You may want to look at using the shape parameter and see if that makes a difference.

Below are four plots of different variants of yours:

Top left: Yours.
Top right: a variation that uses slightly different x values.
Bottom left: Yours at twice the frequency.
Bottom right: a variation that uses slightly different x values at twice the frequency.

As you can see, in all four examples the red and the black overlay exactly.

Yes, in general, these will not be the same. Below is an example with a sawtooth wave, and the two do not line up.

Here is the specific example you ask about with the triangle wave.

R Code Only Below

#43794
#N = 100;
N <- 100
#dsfactor = 3;
dsfactor <- 3
#H = dsfactor*3;
H <- dsfactor * 3
#x = linspace(0,10,N);
x1 <- seq(0,10,10/(N-1))
x2 <- seq(0,10*(N-1)/N,10/N)
#y = sind(18*x);
y1 <- sin(18*x1/180*pi)
y2 <- sin(18*x2/180*pi)
y3 <- sin(36*x1/180*pi)
y4 <- sin(36*x2/180*pi)
#h = ones(H,1);
h <- rep(1,H)

#convyh = conv(y,h);
convy1h <- filter(y1,h)
convy2h <- filter(y2,h)
convy3h <- filter(y3,h)
convy4h <- filter(y4,h)

downsample <- function(signal, factor)
{
  return(signal[seq(1,length(signal), factor)])
}

#convyh_ds = downsample(convyh,dsfactor)./H;
convy1h_ds <- downsample(convy1h, dsfactor)/H
convy2h_ds <- downsample(convy2h, dsfactor)/H
convy3h_ds <- downsample(convy3h, dsfactor)/H
convy4h_ds <- downsample(convy4h, dsfactor)/H
#convyh_ds3 = conv(downsample(y,dsfactor),downsample(h,dsfactor))./(ceil(H/dsfactor));
convy1h_ds3 <- filter(downsample(y1,dsfactor), downsample(h,dsfactor))/ceiling(H/dsfactor)
convy2h_ds3 <- filter(downsample(y2,dsfactor), downsample(h,dsfactor))/ceiling(H/dsfactor)
convy3h_ds3 <- filter(downsample(y3,dsfactor), downsample(h,dsfactor))/ceiling(H/dsfactor)
convy4h_ds3 <- filter(downsample(y4,dsfactor), downsample(h,dsfactor))/ceiling(H/dsfactor)
#figure;
#plot(convyh_ds,'DisplayName','downsample(conv)','LineWidth',1);
#hold on;
#plot(convyh_ds3,'DisplayName','conv(downsample)','LineWidth',1);
#grid minor
#legend('show');
#line([0 length(convyh_ds2)],[0 0],'LineStyle','--');

par(mfrow=c(2,2))
plot(convy1h_ds)
lines(convy1h_ds3, col='red')

plot(convy2h_ds)
lines(convy2h_ds3, col='red')

plot(convy3h_ds)
lines(convy3h_ds3, col='red')

plot(convy4h_ds)
lines(convy4h_ds3, col='red')


How to find tennis courts in aerial imagery I'm interested in finding all tennis courts (and other similar well defined features like basketball courts) in my county, and I have aerial imagery of good (but varying) resolution, but I'm not sure of the best way to find them. Here are two examples of the imagery:

I've looked at the various methods, and I think template matching wouldn't work as it would be very slow since there can be arbitrary scale and rotation, and also the color can vary. The Hough transform sounds promising, but once I get all the lines I'm not sure how to find lines that constitute a rectangle with the appropriate ratio (about 36x29 feet), or better yet to account for the other marked lines.

For background, I'm aiming to add all tennis courts in my county to OpenStreetMap.

You have some very strong color and geometry cues you can leverage. I would try the following:

Extract the Green channel & apply watershed type algorithm on it, followed by connected components. Subsequently compute component statistics (area & bounding box) for each component. Retain only the components with area ~= bounding box size. This will be true only for rectangular objects and will eliminate forests/wooded areas etc.
Isolate the white channel (R=G=B) and apply hough transform on the output. This will give you the lines.
Combine 1 & 2 to get your tennis courts.


Polar form of the Fourier transform of $\sin(t)$ I came across the Fourier transform of sin(t).

It ends up being a purely imaginary (dirac delta) impulse pair. But when considering the frequency domain representation of a signal, we consider the graphs of (i) the magnitude of the signal w.r.t. frequency and (ii) the phase (i.e. argument or angle) of the signal w.r.t. frequency.

In the case of $\sin(t)$, the phase can be seen to be 90 degrees, treating the delta functions in manner analogous to real numbers.

But how would I find the magnitude? Could someone please provide a reference to a textbook?

As another example, what would be the magnitude and phase of

$$\operatorname{FourierTransform}(\cos(t)) \cdot (4 + 2j) \text{ ?}$$

I read two question here. One is about the graphical representation of the Fourier transform of a sinusoidal function, and the other is about the magnitude (and phase) of such a Fourier transform. For ordinary functions, these two would be equivalent, but for Fourier transforms containing Dirac delta impulses they are not. The reason is that mathematically you can't compute the magnitude of a Dirac delta impulse. It is not a function, so there are no function values to take the magnitude of.

A standard graphical representation of the Fourier transform of the function $\sin(\omega_0t)$ would usually just be two arrows at $\omega=\omega_0$ and $\omega=-\omega_0$, one pointing up, the other one pointing down, where the $y$-axis is labeled by $j$ (which makes a separate phase plot unnecessary). An example of such a graph is shown here.

Even though you can't compute the magnitude of a Dirac impulse, you could write down a representation of the Fourier transform of a sinusoid consisting of a sum of Dirac impulses with positive areas and a phase term. For the function $f(t)=\sin(\omega_0t)$ such a representation would be

$$F(\omega)=\pi[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]e^{-j\frac{\pi}{2}\text{sign}(\omega)}\tag{1}$$

(Note that I use the non-unitary definition of the Fourier transform, which is different from the unitary version used by Wolfram Alpha. This is why there's a factor of $\pi$ in $(1)$ instead of $\sqrt{\frac{\pi}{2}}$.)

A similar representation of your example $A\cdot\mathcal{F}\{\cos(\omega_0t)\}$ with complex $A=|A|e^{j\phi_A}$, where $\phi_A$ is the argument of the complex number $A$, is given by

$$A\cdot\mathcal{F}\{\cos(\omega_0t)\}=|A|\pi [\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]e^{j\phi_A}\tag{2}$$


Amplitude demodulation of sinusidal signal in matlab AM-demodulation task and I was asked to demodulate the signal in MATLAB and the result should contain three pictures.

plotting the $fft$ the (frequency spectrum)of the signal i was able to determine that the pictures are located at 0.1 , 0.25 and 0.5 frequencies.

The problem is that my code below only plot one picture which is located at 0 and I don't know how plot the shifted pictures that are located at 0.25 and 0.5 frequencies.

I tried to multiply it by complex exponential and still not able to plot the pictures

The main problem is that i need to shift those pictures that are located at 0.25 and 0.5 frequency to the center and apply my low pass filter, however i tried as the following code but failed. please help me.

The is the $fft$ of the signal.

  I get a mysterious signal modulated with the message as a signal vector z[n],
  by running the MATLAB function mkhwdata. The parameters of this function are your
  social security number in one
  column vector, e.g.
  z = mkhwdata ([xxxxxx]);
  The response from the function will be that a vector z containing the signal z[n], where the vectors
  first element corresponds to the time index n = 0.

  • To convert a vector, e.g. vector s describing the signal s[n], to an image use the provided function present_image. Try for example to watch the signal z[n]
  directly:
  present_image (z)

clc
close all;
clear

fc = 0.2;
fs = 0.5;
[z]= mkhwdata([xxxxxxxxx]);

[Z,nu] = tdftfast(z);
%N = length(z);
 n = [-0.5:1/0.5:0.5];
%n = N*(0:N-1);
 mult  = 0.5.*exp(1i*2*0.25*pi.*n)  + 0.5.*exp(-1i*2*0.25*pi.*n);
 mult2 = 0.5.*exp(1i*2*pi*0.5.*n)   + 0.5.*exp(-1i*2*pi*0.5.*n);
 [z1] = ifft(Z.*mult);

 [z2] = ifft(Z.*mult2);

 %Lowpass filter
 [num,den] = butter(9,0.2);

 y  = filter(num,den,z);
 y1 = filter(num,den,z1);
 y2 = filter(num,den,z2);
 figure(4)
 ylabel('Amplitude');
 xlabel('Frequency (hz)');
 subplot(221)
 present_image(z);
 subplot(222)
 present_image(y);
 subplot(223)
 present_image(y1);
 subplot(224)
 present_image(y2);

I think there are some problems with your analysis structure.

First of all, for being able to see from -0.5 to 0.5, without aliasing, you have to sample at 1.

Secondly, if you want to down convert a signal, bring it to DC you need to multiply it with the corresponding cosine in time, or discrete time, or apply the convolution to its spectrum with the corresponding Fourier transformed cosine (a cosine in spectrum is a delta Dirac function on the frequency of the same cosine).

By the way, I do not know what the mkhwdata function does.
The solution should look something like this:

clc
close all;
clear

fc = 0.2;
fs = 1; %%%%%%%%%%%%%%%%%%

[z]= mkhwdata([xxxxxxxxx]);

[Z,nu] = tdftfast(z);
%N = length(z);
n = [0:(length(z) - 1)];
%n = N*(0:N-1);
mult  = exp(1i*2*0.25*pi.*n);
mult2 = exp(1i*2*pi*0.5.*n);

% We bring the image to DC
[z1] = (z.*transpose(mult));
[z2] = (Z.*transpose(mult2));

%Lowpass filter
[num,den] = butter(9,0.2);

y  = filter(num,den,z);
y1 = filter(num,den,z1);
y2 = filter(num,den,z2);

figure(4)
ylabel('Amplitude');
xlabel('Frequency (hz)');
subplot(221)
present_image(z);
subplot(222)
present_image(y);
subplot(223)
present_image(y1);
subplot(224)
present_image(y2);


Is there better than hamming code for correcting errors As known, hamming code can be represented by (7,4) where data are 4 bits, and 3 bits for parity,  and (15,11) where 11 bits are data and 4 bits for parity. It can correct up to one bit, or detect two bits without correcting.

Parity bits will cause an issue of spectral efficiency or in data rate, so my question, is there another code better than hamming code which can correct more errors with less parity bits ?

I think the question hinges around the trade off of how heavily you code the information and the rate. Generally adding more parity bits will give better error detection/correction ability but it will reduce the communication rate since each each information word is being mapped into a longer code word. Because of this trade off, I am not sure you will find better codes than the mentioned code in terms of it having less parity bits while increasing error detection/correction ability.


How do I get the Time information after STFT? I'm tring to code a STFT program in Matlab. I got results but I'm confused a bit.
I have a audio signal of 900000 pts sampled at 1000Hz (15 min signal).
I've taken fft with following parameters =>

hop size of H = 1400 pts
the window length M = H * 6

So, 5/6 of the window overlaps with the next one.
The program runs for 637 iterations, as (637*1400 + M = 900200)

Now, I get a 637 pt vector with frequency in them. I should get a frequency vs time plot. But it only has 637 frequency bins, for a signal of 900000 time sample.

How, can I find the change of frequency with Time (900000 sample) or
15 min?
Moreover, how does changing the window length and overlap affect the outcome of STFT? The code is below.

%hop size
L = 1400;

%window length & Formation
D = 6;                        %User Defined Parameter in Cooper's paper
M = L*D;
w = hanning(M, 'periodic');    %window formation

k = 1; % k =index
j = 1; % j =result vector index
while k + M <= xlen

  %data chunk
  %x = x(k:1400);

  %window length & Window
  xw = x(k+1:k+M) .*w;                  %windowed signal

  %zero padding
  b = 5;                           %Zero Padding Factor 'b'
  xw = [xw; zeros(M*5,1)];
  [nft,c] = size(xw);

  %plot(abs(fft(xw)))
  X = abs(fft(xw));
  XF = X(1:length(X)/2);
  [q,ind] = max(XF);
  v(j) = ((fs/2)/(nft/2)) * (ind-1);  %convert Bin to frequency
  k = k+L;
  j = j+1;
end

figure(1);
plot(v)
axis([0 1400 0 100])

figure(2);   %zoom version
plot(v)

How, can I find the change of frequency with Time (900000 sample) or 15 min?

OK, so you need 900000 frequencies.  Well, one way to do that is to take each element of your 637 length vector and expand it to (up to) $H$ samples.

  Moreover, how does changing the window length and overlap affect the outcome of STFT?

The longer the window length, the better the frequency resolution is. The shorter the length, the worse the frequency resolution will be.

Before we can answer your question, we will need to know what form your signal takes.

If it's a single tone in noise, then you may be able to use one of these or one of these or this approach.


Relationship between z-transform and DFT I'm studying for a Signals Processing exam and came across an exercise that I'm finding pretty difficult to solve. It says:

  Asume there is a signal $x[n]$ of length $N$. Its $\mathcal{Z}$-Transform is $X(z)$. Find the expressions of the DFTs of the next signals using a
  correct sampling of $X(z)$.

  a) $x_{1}[n]=x[n]$  if  $0\le n \le N-1$ and $x_{1}[n]=0$  if  $N\le n \le 2N-1$

  b) $x_{1}[n]=x[n]+x[n-N]$

For a), I thought of expressing  $X[k] = X(z) \bigg|_{z=e^{j\frac{2\pi}{N} k}}$, i.e. the $\mathcal{Z}$-Transform evaluated in $z=j\frac{2\pi}{N} k$ is the DFT of $x[n]$. Due to the zeros added at the end of the original signal, I arrived at the conclusion that  $X[k] = X(z) \bigg|_{z=e^{j\frac{2\pi}{2N} k}}$. Namely, the $\mathcal{Z}$-Transform is now sampled in $2N$ points around the unity circle. Is this reasoning right?

And regarding exercise b), I have absolutely no idea how to solve it.

Thanks for your time!

For a) you're correct. For b), $x_1$ is a length $2N$ signal, and its DFT is given by

$$X_1[k]=\sum_{n=0}^{2N-1}x_1[n]e^{-j2\pi kn/2N}=\sum_{n=0}^{2N-1}x_1[n]e^{-j\pi kn/N}\tag{1}$$

With $x_1=x[n]+x[n-N]$ you get

$$\begin{align}X_1[k]&=\sum_{n=0}^{N-1}x[n]e^{-j\pi kn/N}+\sum_{n=N}^{2N-1}x[n-N]e^{-j\pi kn/N}\\&=X\left(e^{j\pi k/N}\right)+\sum_{n=0}^{N-1}x[n]e^{-jk\pi(n+N)/N}\\&=X\left(e^{j\pi k/N}\right)+e^{-jk\pi}X\left(e^{j\pi k/N}\right)\\&=X\left(e^{j\pi k/N}\right)(1+(-1)^k)=\begin{cases}2X\left(e^{j\pi k/N}\right),&k \text{ even}\\
0&k\text{ odd}\end{cases}\tag{2}\end{align}$$

Eq. $(2)$ shows that decimation in the frequency domain corresponds to aliasing in the time domain.


How is multiplexing on FM signals performed? How does multiplexing differ from FM signals to AM signals. If both FM and AM shift the message signal in frequency, then  is there a difference between multiplexing for FM or AM signal?

In U.S. broadcast FM, the various subcarriers (pilot tone, stereo L-R, RDS, etc.) are modulated into the upper sidebands (above the 15 kHz L+R/mono audio, but below 100 kHz, for instance 38 kHz DSB for the stereo subcarrier) of the baseband signal before that entire multiplex is used to wideband frequency  modulate a single, much higher frequency, FM carrier.

The FM carrier is modulated, but not shifted in center frequency.


Plotting Spectrum Analysis I am working on a spectrum analyzing simulator, which is meant to compute and plot the FFT of a sine wave carrier signal of $16 Hz$.

From my practicals with a spectrum analyzer, I know that the spike shows up on the $16Hz$ mark. But after computing my own FFT and passing the magnitude of the FFT to a plotter in Matlab, I get two spikes, one at the beginning and the other at the end.

My sample rate is $2048$ samples per second, and what I plotted was the respective $abs(FFT)$ against a frequency step of $16/2048$. Doing the fftshift only moves the spike to the $8Hz$ mark and not $16Hz$.

I've tried computing FFT in different programming languages and with different tools (java jtransform, Matlab, C language, Excel, two java written classes for FFT by different universities).

I dont think my FFT result are wrong, I think it has to do with the plotting. I do not understand why my plot does not look like it would in real life. From my description and information, can you see any problems that would cause my plot to be wrong?

Generic FFTs use complex-valued input arrays and provide complex-valued output arrays, so they can differentiate positive and negative frequencies for complex-valued signals.   For real-valued signals, which is what it sounds like your case may be, this means that the positive and negative frequencies are mirrors of each other, so the spike at 16Hz will be mirrored by a spike at -16Hz.  Matlab's FFT function output array has the positive frequencies first followed by the negative frequencies, so one would expect two spikes:  the first at 16Hz, and the second wrapped to -16Hz which will be almost at the end of the array.

fftshift should have "unfolded" the output array so that the negative frequencies plot on the left, with DC (0Hz) in the center, and positive frequencies on the right.   You should still see two spikes, one at -16 Hz and one at +16Hz, symmetric about zero.


The description of the MRI signal from the magnetization vector I have been working on simulating the MRI signal and what I have done so far is simulate the actions of the RF pulse and the required localization gradients.
So, in the simulation, I am looking at a gradient echo sequence and I have the RF pulse and the slice selection gradient which tips the spin magnetization in the XY plane. This is followed by a rephasing gradient (to take into account the dephasing due to the slice selection gradient).

After this, I have the PE gradient along the Y-direction and the frequency rewinder (to have a symmetrical echo) and the readout gradients along the X-direction.

In the simulator, I am tracking the magnetization vector during the whole process by solving the Bloch equations. Now, the part that I am confused about is how can I convert this to the MRI signal? At the end of the day for each spin, I have its magnetization at every time step. During the readout process, I would like to convert it to the MRI signal and I am not sure how to go about it.

You always measure the (vector) sum over all magnetization vectors in your object, $S(t) = \sum\limits_{(x,y)} M(x,y,t)$. This sum gives you an amplitude and a phase. That is your measured signal.

From the gradients that you simulate, you know the $k$-space position as a function of time, i.e. you have $k(t)$. Now you sample your signal at certain points and sort the signals $S(t_n)$ into the respective $k(t_n)$-space. A FFT should now convert your $k$-space into the image again.


What is the advantage of applying a Hamming window for MFCC calculations? I am trying to implement MFCC by following this and I found that there is a step in which signals are divided into frames and then passed on to Hamming windowing process. The reason given there is to correct the discontinuity at the start and last samples of the frame. OK that sounds good, then the output of the window is fed for FFT process. So I thought of doing this sample code separately.

My code:-

%Creating a signal to perform hamming windowing on.
fs=8000;
ts=1/fs;
N=1000;
t=(0:ts:(N-1)*ts);
x=sin(2*pi*60*t);

figure,subplot(221),plot(t,x,'.-'),title('Original Speech'),xlabel('Time');grid on
f=(-fs/2:fs/(N-1):fs/2);
subplot(223),plot(f,fftshift(abs(fft(x)))),title('Frequency Spectrum'),xlabel('Freq (Hz)');grid on
windowed=x.*hamming(length(x))';
subplot(222),plot(t,windowed,'.-'),title('Hamming Window Applied'),xlabel('Time');grid on
subplot(224),plot(f,fftshift(abs(fft(windowed)))),title('Frequency Spectrum After Hamming'),xlabel('Freq (Hz)');grid on

The plot generated is:-

My Problem is:-

What is the advantage of Hamming Windowing. I am a noob, so according to me, firstly we are loosing the amplitude at the staring and ending of the signal in the time domain and secondly, in the frequency domain, there is very little difference. So what is the advantage of doing all this much computation.
Please explain with as much detail as possible.

And it would be nice if someone could point out to more detailed (according to them) tutorial then what I am referencing or what they read while learning MFCC themselves.

Ideally you would only see that frequency as a narrow peak but because of the finite-length window you also get that other artificial crud: Windowing equals time domain multiplication by the window function (here rectangular or Hamming). Time domain multiplication equals frequency domain convolution: All frequencies will be replaced by the Fourier transform of the window function. This is what those look like:

The convolution spreads the true frequency to frequency bins around it. By using a window function other than the rectangular one there is less of the sidelobes compared to the main lobe, so less artificial long-distance spread making the result cleaner and better suited for frequency-selective analysis. You seem to have by accident chosen a frequency where you don't see the sidelobes with the rectangular window because the dips between the sidelobes land between your frequency bins. So you just see slopes.


Band limited wavetables I am a bit confused with the concept of band limited wavetables and the methods such as BLET, BLIP, MinBLEP, polyBLEP etc.

As far as I understand certain waveforms contain infinite harmonics, which alias back into the signal, even if they were generated at a certain frequency below the Nyquist rate. According to this paper Alias-Free Digital Synthesis of Classic Analog Waveforms, these basic waveforms have aliasing "due to having to round off the discontinuity time to the nearest available sampling instant".

What and were is the discontinuity of the waveform if the wavetable contains e.g. a full cycle of a sawtooth wave at a certain frequency? (where looping the table will restart the new cycle from the last point of the old one)

After reading through some code on how MinBLEPs work, Im confused with how all these steps are linked. E.g. the signal is windowed and FFT'd and I assume the cepstrum function is used to look for the discontinuity inside the FFT bins. I dont get what the minimum phase is used for. The signal is then run through the iFFT.

I please need help in understanding how all these techniques are combined and what exactly is done in the practical sense.

After all of this, is the minBLEP a kind of wavelet that is stored in a buffer and then applied to the signal being generated, so that it smooths the discontinuities? At first I thought this was all done inside an oscillator which needs some level of look ahead to constantly detect discontinuities and apply the minBLEP wavelet at the correct time, but now I don't get how this is done inside the wavetable generator function itself (offline), that an oscillator can just scan in real time.

It looks like you are mixing up two techniques: MinBLEP and band-limited wavetables.

Band-limited wavetables

Let's assume you have a 256 sample long array with a waveform of interest - for example a naive sawtooth, with the first sample in your array at -1 and the last sample at 1.

To replay this data at a variable pitch, two techniques are possible:

First, you can control the rate at which the data is played. For example, to get a 440 Hz pitch, you replay the 256 samples verbatim, at a rate of 112640 Hz. This is possible with custom hardware, and is the technique used in early wavetable synthesizers (PPG Wave, Prophet VS). This requires one clock/DAC per oscillator, and is rather expensive!

Or you use a fixed sample-rate, and you modify the pitch by interpolation through the wavetable. Let's say the sample rate is 48kHz. Playing the data one sample after the other yields a tone with a frequency of $48000 / 256 = 187.5$ Hz. To play a tone at 440Hz, you would need to play the data at $440 / (48000 / 256)=2.35$ its original rate, and that's where you get into trouble. Two ways of understanding why:

You can consider the discontinuity. You'll reach the end of the table and "wrap around" between the 109th and 110th sample. So the discontinuity should actually be located between these samples.
Or you can think in terms of frequency content: by playing the data at 2.35 its original rate, any harmonic it contains at a frequency above nyquist / 2.35 will be aliased.

There are two ways of dealing with the problem. A first approach is to use band-limited interpolation, which is a way of replaying the data in such a way that aliased components are filtered. Another solution is to pre-compute several variants of your waveform in which the highest harmonics are progressively filtered. For example, when you have to playback the data at twice its original rate, you use a different sawtooth wavetable in which all harmonics above nyquist/4 are filtered-out. Obviously, you can't have a table ready for every possible playback rate, so a common solution is to pre-compute a dozen of such tables (for example one per octave), and interpolate between them. Korg did that on their DW series, I used it on my DIY AVR-based synths.

There are different ways of building such a band-limited table. You can work by additive synthesis (adding the right number of harmonics), or try to derive the band-limited waveforms by operations on a band-limited impulse (sinc waveform).

minBlep

This technique is completely unrelated to wavetable synthesis. The idea is to precompute a table containing a band-limited and causal step waveform. This computation is done once and for all, and involves an FFT/IFFT - that's probably where you encountered these operations, which are used only at this offline step. The code you link to in your question only covers this step!

Then, in real-time, you generate the naive waveform, and whenever there is a discontinuity in the waveform, you add, over the next 40 samples or so, samples from the band-limited step table to smooth the edge. Which samples you need to add from the minBlep table are determined by the fractional sample at which the discontinuity occurs. No look-ahead is involved: the point of minBLEP is that the correction is causal, so you don't have to start adding it before the discontinuity occurs.

Note that the computational cost of this method will increase as the frequency goes up, since the discontinuities will be more frequent - you'll have to keep track of an increasing number of overlapping corrections.

I wouldn't recommend minBLEP by the way. PolyBLEP uses a polynomial approximation of the band-limited step function and yields very good results. No lookup data is required and the code is easier to write because the correction only lasts for a handful of samples.

The two together?

It would be a strange idea to mix the two techniques together, but someone could use minBLEP offline to precompute a band-limited sawtooth for various band-limiting ratios. For example, a 256 sample-long naive sawtooth, then a 256 sample-long sawtooth with no harmonic above 0.25, then a 256 sample-long sawtooth with no harmonic above 0.125... And then one could use wavetable playback to replay them - selecting/crossfading one of these tables depending on the required playback rate. This would be a bit odd, there are simpler ways of building wavetables!


Recovering an image using IDCT through DCT coefficients I have computed the discrete cosine transform (DCT) of an image in MATLAB. I got a matrix of discrete cosine transform coefficients in the workspace. How can I recover the original image from its discrete cosine transform coefficients in MATLAB?

If you used dct2  on an image, or on each image channel separately for color images,  with the appropriate input format, then idct2 recover the original image (the transformation is invertible and orthogonal), up to rounding errors (similar discussion in Why is $\exp(\ln(x))−x\neq 0$ in floating point arithmetic?).

You can find for instance a demo below:

dataRGB = double(imread('autumn.tif'));
data = rgb2gray(dataRGB/255);
dataDct=  dct2(data);
dataDctIdct = idct2(dataDct);

subplot(2,2,1)
imagesc(data);
xlabel('Original (grayscale)');
subplot(2,2,2)
imagesc(log(abs(dataDct)));
xlabel('Original DCT');
subplot(2,2,3)
imagesc(dataDctIdct);
xlabel('Original recovered');
subplot(2,2,4)
imagesc(data-dataDctIdct);
xlabel('Difference');
colorbar
colormap gray


Composing two Impulse Responses Pardon my layperson understanding. My application is audio design.

What is the correct operation / algorithm to compose two impulse responses (IRa, IRb) into a third impulse response IRc such that it equals the response of running a signal through IRa then IRb?

Am I correct that since IRs represent the response to the delta function, that I can simply convolve IRa by IRb directly? i.e:

IRc[n] = IRa[n] * IRb[n]

Without needing to worry about deconvolving the output of the first IR convolution?

Thank you.

Yes, that's correct.

The impulse response of a cascade of two systems is the convolution of the individual Impulses responses. Convolution is commutative, so it doesn't matter in which order you convolve.

Impulse responses are ONLY defined for linear time invariant systems, so this will work for systems like equalizers or filters but NOT for things like dynamic processors or limiters.


one to many images matching in matlab? i have database for 50 person ,it is 250 images (5 images/person) , i need to perform one to many matching between images for recognition persons.

so , i divided the image as 2 parts (one for training and one for testing)

for example if we say that the first folder contains 100 images(2 images /person),and the second contains 150 images (3 images /person)

is this division correct?
 if it is correct , how can i matching each image in folder1 with all image in folder 2(which will contain 3 similar image) and obtain the better result after (that give max matching percentage)from folder 2 and finally returning the name of this image ?

i need to perform this in matlab

please , i need many advices and consultations with any one has good knowledge in this field.

this is an example of my images:

regards

5 images per person is really very few.  In this situation it makes sense to use as much data for the training set as you possibly can.  This means for each person randomly select 4 images for training, and one for testing.  You should use this procedure several times so that you have a number of training set/test set combinations.

If you have access to the Statistics Toolbox for Matlab, check out cvpartition function.

If not, then you should have a loop over all persons, and for each person use the randperm function to generate a randomly shuffled list of indices of the images.  Then you can take the image corresponding to the first index and add it to your test set, and put the rest into your training set. Or you can put 2 images into the test set, or however many you like.  Generally, the training set should be larger than the test set.

You should generate the training and the test set several times, each time calling rng with a different value, to reset the seed of the random number generator.  If you generate 10 training/test set pairs, then you evaluate your algorithm 10 times and get 10 accuracy values.  From them you can compute error bars, or you can compare to other algorithms using hypothesis testing.

As to how you actually classify the images, that is very much a separate question, and the answer depends on what's in the images.  What exactly is on that image you've posted?


Inverse $\mathcal Z$-transform of rational functions What will be inverse $\mathcal Z$-transform for this function:

$$H(z) = \frac{\left(1+\beta z^{-1}\right)\left(1+\beta z\right)}{\left(1+\alpha z^{-1}\right)\left(1+\alpha z\right)}$$

Given the $\mathcal Z$-transform :
$$H(z) = \frac{\left(1+\beta z^{-1}\right)\left(1+\beta z\right)}{\left(1+\alpha z^{-1}\right)\left(1+\alpha z\right)}$$

The inverse $\mathcal Z$-transform can be found by the method of partial fraction expansion guided with the Region of Convergence ROC associated with $H(z)$

I will proceed with the most fundamental method and apply necesary algebraic manipulations to get the required standard form as follows:

\begin{align}
H(z) &= \frac {\left(1 + \beta z^{-1}\right)\left(1 + \beta z\right)}{\left(1 + \alpha z^{-1}\right)\left(1 + \alpha z\right)} &\scriptstyle{\text{convert to negative power of $z$}}\\
 &= \frac{\left(1 + \beta z^{-1}\right)\left(z^{-1} + \beta\right)}{\left(1 + \alpha z^{-1}\right)\left(z^{-1} + \alpha\right)} &\scriptstyle{\text{convert to standard form of } (1-d_k z^{-1})}\\
&= \frac{\beta}{\alpha} \frac{\left(1 + \beta z^{-1}\right)\left(1 + \frac{1}{\beta} z^{-1}\right)}{\left(1 + \alpha z^{-1}\right)\left(1 + \frac{1}{\alpha}z^{-1}\right)} &\scriptstyle{\text{expand the brackets}}\\
 &= \frac{\beta}{\alpha} \frac{\left(1 + \left(\beta + \frac{1}{\beta}\right) z^{-1} + z^{-2}\right) }{\left(1 + \left(\alpha + \frac{1}{\alpha}\right) z^{-1} + z^{-2}\right)} &\scriptstyle{\text{apply the long division}}\\
&= \frac{\beta}{\alpha} \left( 1 + \frac{ \left(\beta + \frac{1}{\beta} - \alpha - \frac{1}{\alpha}\right)z^{-1} }{\left(1 + \alpha z^{-1}\right)\left(1 + \frac{1}{\alpha}z^{-1}\right)} \right) &\scriptstyle{\text{go further...}}\\
&= \frac{\beta}{\alpha} \left( 1 + \frac{A}{1 + \alpha z^{-1}} + \frac{B}{1 + \frac{1}{\alpha}z^{-1}} \right) &\scriptstyle{\text{where $A$ and $B$ are found from the right side quotient of $H(z)$}}\\
\end{align}

Now assuming $A$ and $B$ are found as functions of $\alpha$ and $\beta$, we shall define all possible ROC of $H(z)$ to find the corresponding inverse transform $h[n]$.

Now If $0<|\alpha|<1<\frac{1}{|\alpha|}$ and ROC is such that $|z|<|\alpha|$ then: $h[n]$ will be unstable and anti-causal as
$$h[n] = \frac{\beta}{\alpha} \delta[n] + \frac{\beta}{\alpha} \left(- A\left(-\alpha\right)^n \text{u}[-n-1] - B\left(- \frac{1}{\alpha}\right)^n \text{u}[-n-1] \right)$$
Or else if $0<|\alpha|<1<\frac{1}{|\alpha|}$ and ROC is such that $|\alpha| <|z|<\frac{1}{|\alpha|}$ then: $h[n]$ will be stable and two-sided as
$$h[n] = \frac{\beta}{\alpha} \delta[n] + \frac{\beta}{\alpha} \left( A(-\alpha)^n \text{u}[n] - B\left(- \frac{1}{\alpha}\right)^n \text{u}[-n-1] \right)$$
Or else if $0<|\alpha|<1<\frac{1}{|\alpha|}$ and ROC is such that $|z|>\frac{1}{|\alpha|}$ then: $h[n]$ will be unstable and causal as
$$h[n] = \frac{\beta}{\alpha} \delta[n] + \frac{\beta}{\alpha} \left( A(-\alpha)^n \text{u}[n] + B\left(- \frac{1}{\alpha}\right)^n \text{u}[n] \right)$$
The cases for $0<\frac 1\alpha<1<\alpha$ are simply the same, merely interchanging $\alpha$ with $\frac 1\alpha$