Why complex numbers are represents as a+ib and can't be as (a,b)? I am confused as why do we need to represent the complex numbers with the imaginary y-axis if we can simply represent them as (x,y) ?

I've read that Multiplication by i is an anti-clockwise rotation of a quarter-circle over y-axis.

Multiplying 1 by i gives i. Multiplying i, by i once more, does another quarter-circle and gives -1. So, multiplying by -1 means a rotation of a half-circle. That is the meaning of i*i=-1.

So what does that supposed to mean ?

Suppose, I am solving an equation and I ended up with an answer like 3i does that mean i've moved from x-axis to y-axis by half-circle counter clock wise ? I couldn't be able to visualize this properly

Yes, in signal processing, complex numbers are usually visualized on the complex plane, as you have said.

The reason is that if you put them on a plane, then you are able to measure two important quantities:

1) Magnitude, which is $\sqrt{x^2 + y^2}$

2) Phase angle between your point and the origin, given by $\tan^{-1} \frac{y}{x}$.

If you simply left them as a point, ($x$,$y$), you would not be able to concretize and have a frame work for those quantities.

You may ask, why are those quantities, in turn, important? In signal processing, we are of course dealing with signals, and physically, we are dealing with 'real' signals. However, though a nice trick, an constant oscillation of a quantity in 'real' life, (like a cosine wave), is equivalent to two phasors, rotating around in opposite directions on the complex plane, and adding up together. With this framework, we can see that the phase angles 'cancel' each other out, and that the magnitudes of their resultant give us the magnitude of our 'real' signal.

In fact this is what one of euler's formulas captures. That is:

$$
\cos(2\pi ft) = \frac{e^{j2\pi ft} + e^{-j2\pi ft}}{2}
$$

You can see here how we can easily relate a 'real' world concept, like an oscillating cosine wave, with the 'complex' world of phasors, as they exist and rotate around in the complex plane.

This is one of the corner stones of DSP.


Confusion for two-sided signal Given a continuous LTI system with transfer function
$$H(s)= -\frac{2s}{(s+6)(s+2)}$$

Plot the location of the pole(s) and zero(s)
Find all possible regions of convergence
From the problem above find the impulse response

Here's my attempt

The poles are $s_1=-6$ and $s_2=-2$ and the zero is $s_3=0$

The all possible roc-s are

Using inverse laplace transform I found $h(t)=e^{-2t}-3e^{-6t}$. For the right-sided signal I multiply $h(t)$ with $u(t)$ and I got $h_{RS}(t)=(e^{-2t}-3e^{-6t})u(t)$ and for the left-sided signal I multiply $h(t)$ with $u(-t)$ and I got $h_{LS}(t)=(e^{-2t}-3e^{-6t})u(-t)$.

My confusion is at the two-sided signal. My solution is $h_{TS}(t)=e^{-2t}-3e^{-6t}$ and my friend's solution is $h_{TS}(t)=-e^{-2t}u(-t) - 3e^{-6t}u(t)$. When I plotted the signal both are a two-sided signal. Since my professor didn't submit the solution we didn't know who's right or wrong.
Thank you.

Given LTI system with transfer function :
$$H(s) = \frac{ -2 s}{(s+6)(s+2) } \tag{1}$$
which has the poles $s = \{ -2 , -6 \} $, and zeros $ s = 0 $.
Apply partial fraction expansion :
$$H(s) = \frac{1}{s+2} - \frac{3}{s+6} \tag{2}$$
Based on pole locations, there will be three ROCs,

ROC1 :  -$2 < \mathcal{Re}\{s\} < \infty $  ; stable
ROC2 : $-\infty < \mathcal{Re}\{s\} < -6$  ; unstable
ROC3 : $ -6 < \mathcal{Re}\{s\} < -2$   ; unstable

The impulse responses found from (2) and the chosen ROC , with the fact that
$$\mathcal{L}\{ e^{at}u(t) \} \longleftrightarrow \frac{1}{s-a} ~~,~~ \mathcal{Re}\{s\} >\mathcal{Re}\{a\} \tag{3}$$
$$\mathcal{L}\{ -e^{at}u(-t) \} \longleftrightarrow \frac{1}{s-a} ~~,~~ \mathcal{Re}\{s\} <\mathcal{Re}\{a\} \tag{4}\\\\$$

ROC1 :  -$2 < \mathcal{Re}\{s\} < \infty $  ;  right-sided  : $h(t) = (e^{-2t}-3 e^{-6t}) u(t)\\\\$
ROC2 : $-\infty < \mathcal{Re}\{s\} < -6$  ; left-sided : $h(t) = (-e^{-2t}+3 e^{-6t}) u(-t)\\\\$
ROC3 : $ -6 < \mathcal{Re}\{s\} < -2$   ; two-sided : $h(t) = -e^{-2t}u(-t)-3 e^{-6t}u(t)$.


Save computational burden by calculating the output of decimate filter at output rate? I read in a book that the output of a decimate filter can be calculated at output sampling rate. In this way the number of filtering operations can be greatly reduced. But I am wondering how to implement that in practice.
The built-in filter functions provided by Matlab and a lot of realtime DSP libraries are optimized to have best computation and storage performance, but they don't provide the capability of calculating the filter output only at output sampling rate. Is it possible to write a code, in a for-loop style, that only calculates the decimated data, and outperforms the performance of built-in functions? If not, what modifications can be done to the built-in functions to reduce the computation of decimation filter, especially in realtime implementations?
Thank you

Is it possible to write a code, in a for-loop style, that only calculates the decimated data, and outperforms the performance of built-in functions?

Yes.
(the way you asked this, here would be where the answer is over, but it's not what you meant, sooooo let's continue)
"for loop style" isn't very specific, and any algorithm that implements some form of convolution will iterate over elements, but sure:
For decimation by a factor of $M$:
Decompose your filter in $M$ subfilters. Round-robin the input samples into the individual subfilters. Every $M$ input samples, calculate the output as sum of the outputs of the subfilters.
There's many publications on polyphase filters. harris' multirate signal processing is the go-to reference; but there's really many other references, too, see for example this slide deck from the Binghamton University (I'm not even sure the graphics are from that slide deck's author, my apologies if I'm wrong, but I feel I've seen the graphics before).
I recently added a different answer on the same topic, with a bit more background on the polyphase "Magic" that allows the lower rate.

The built-in filter functions provided by Matlab and a lot of realtime DSP libraries are optimized to have best computation and storage performance, but they don't provide the capability of calculating the filter output only at output sampling rate.

That's not true; any real-time DSP framework I can think of has decimating / resampling filters that do exactly that. (example: GNU Radio's decimating FIR filter function; Liquid DSP's firdecim, csdr fir_decimate, the list goes on: this is really a core functionality of every real-time signal processing toolkit.)
I don't know how efficient matlab's implementation is, but upfirdn with an interpolation ratio of 1 and a decimation > 1 should do exactly that.


Optimization FIR filter I have a doubt about this :
Iterative methods for optimization in FIR filter (on the one hand traditional like window,frequency sampling ,WLS,Remez Exchange...,on the other hand evolutionary methods like SA,PSO,GA....) ,it is used for :

Verify that the filter has linear phase  (symmetric).
or the symmetric is satisfied ,but have another problems in this type of
filters should be reduced or eliminated to be used in wider applications.

Please correct me if I am wrong about that.

I'm not sure if I understand your question correctly, but I'll try to make a few things clear. First of all, note that windowing and frequency sampling are not iterative. On the other hand, the Remez exchange algorithm (used in the ParksMcClellan program) is iterative.

All methods for the design of linear phase FIR filters impose the linear phase constraint in the formulation of the design problem, i.e., the filters are always guaranteed to have a linear phase. The difference between the different methods is either the way the optimal solution is computed, or the design criterion, i.e., the way the error to be minimized is defined. Some common design criteria are the least squared error criterion (used in the Matlab function firls.m), or the maximum (Chebyshev) error criterion, where the maximum error is minimized (used in the Matlab function firpm.m).


What are the mathematical tools that I can employ for automatic music composition? I would like to do research on automatic music composition. At the lowest level , I would like to start from discrete frequency levels as states. Here I would like to know the mathematical tools I can use for this? I have heard of markov chains and chaos theory etc used in this purpose. But one professor recently told that graph theory will be helpful in automated music composition. Is it so? Can someone recommend me the books required for starting such a research.

There is an open source on the similar project you may be interested: http://sourceforge.net/projects/amc/
The program generated chords and bass lines from .wav or .midi input.


What do these values in the freeverb diagram mean? I am implementing a reverb based on the Freeverb implementation detailed here. There is also a more in-depth description of the LBCF in question here, but I am not able to find the information I want at that link.

In the diagram given, they use a lowpass-feedback comb-filter, denoted by LBCF(f,d,N). N is clearly the number of samples. I suspect d means the feedback gain, though I'm not sure if a value of .2 makes sense for that. I have no idea what f means.

So here are my questions:

What does d mean?
What does f mean?
What cutoff frequency is used for the lowpass filters?

It's a Lowpass-Feedback Comb Filter, so let's dissect this one at a time:

Comb filter: You have delay line and add the output back to the original signal. That gives you a single echo
Feedback: You take some of the output of the delay line and put it back into the input of the delay line. That gives you a series of echo's that repeat regularly. You need to put a gain < 1 on this, otherwise it will repeat forever. That gain is called the feedback gain and controls how fast (or slow) the echoes decay.
You put a lowpass into the feedback path to model the typical absorption in a room that's higher at high frequencies and lower at low frequencies. You basically make the feedback gain frequency dependent.

$d$ is the pole location of the lowpass filter. That determines the cutoff frequency of your lowpass filter

$f$ is the feedback gain and determines how quickly the repeating echoes decay. That should be a function of the reverb time (at low frequencies) and the delay size


Auto correlation definition My question has to do with the definition of auto correlation/cross-correlation for random processes.
Oppenheim/Schafer (Discrete time Signal Processing, Pg. 815 (Appendix A.2),2nd ed.) define auto correlation in the following way:
$\phi_{xx}[n,m]=E[\mathbf{x_nx^*_m}]=\int_\infty^\infty x_nx^*_m p_{x_n,x_m}(x_n,n,x_m,m)dx_ndx_m$

{$\mathbf{x_n}$} is a random process and $p_{x_n,x_m}(x_n,n,x_m,m)$ is it's joint probability density function.
What is the significance of using the conjugate of $x_m$ in the above formula? Is there a physical interpretation for it?

it's funny, but i have an old copy of O&S and i would have expected them to be consistent and use "$x[n]$" notation instead of "$x_n$" notation for discrete-time signals.

The expression above

$$ \phi_{xx}[n,m] \ = \ E[\mathbf{x_n \ x^*_m}] \ = \ \int_\infty^\infty \int_\infty^\infty \ x_n \ x^*_m \ p_{x_n,x_m}(x_n,n,x_m,m) \ dx_n dx_m $$

is the average value of the pair $\mathbf{x_n \ x^*_m}$ where the average is computed from the joint probability density function $ p_{x_n,x_m}(x_n,n,x_m,m) $

If the process is stationary and ergodic (and i think that any process that is fully stationary is also ergodic) i think this can be simplified a little

$$ \phi_{xx}[n-m] \ = \ E[\mathbf{x_n \ x^*_m}] \ = \ \int_\infty^\infty \int_\infty^\infty \ x_n \ x^*_m \ p_{x_n,x_m}(x_n,x_m) \ dx_n dx_m $$

now the root meaning of "ergodic" is that probabilistic averages can be replaced with time averages and vise-versa:

$$ \langle \mathbf{x_n \ x^*_{n+k}} \rangle = \lim_{L \to \infty} \frac{1}{2L+1} \sum\limits_{n=-L}^L \mathbf{x_n \ x^*_{n+k}} $$

the process is ergodic if all averages, including the one expressed above, are the same, whether they are probabilistic or computed from the data

$$ \langle \mathbf{x_n \ x^*_{n+k}} \rangle = \phi_{xx}[k] = E[\mathbf{x_n \ x^*_{n+k}}] $$


What is the value of the given periodic signal at any time, how to solve it in MATLAB? Firstly hello all, I need help about periodic signals. I have a question as below.

$x(t)$ is a periodic signal and $0.1t^3[u(t) - u(t - 7)]$
describes its one period.
What is the value of $x(t)$ at time $t = 317$?
Hint: Determine one period of the signal and find its values.

How can we solve this in MATLAB ? I tried with the code below, but answer was not correct
x(t) = 0.1 * t^3 * (heaviside(t) - heaviside(t - 7));
ezplot(x);
x(317)

Thank you.

Assuming continuous-time for the probem, for all periodic signals with a period of $T$ , for any integer $m$ , we have :
$$  x( t + mT) = x(t) \tag{1}$$
So, for example with  $T=7$ , and $m=2$ , you can see that $x(15) = x(2\cdot 7 + 1) = x(1)$. The operator that is used to find base argument for a given $t$ value is the modulus :
$$ x(t) = x( \text{mod}(t,T) ) \tag{2}$$
where $\text{mod}(t,T)$ is the modulus of $t$ wrt $T$.
For your example, your period is $T=7$, and you can see that $x(315) = x( \text{mod}(315,7) ) = x(2)$.
Hence the value of $x(t)$ for $t=2$ is $0.1 \cdot 2^3 = 0.8$.
I don't know if there's anything to solve with MATLAB here though?


Use of the harris-Moerder Nyquist Pulse Shaping Filter I became aware today through fred harris' excellent presentation at the DSP Online Conference (https://www.dsponlineconference.com/) of the harris-Moerder pulse shaping filter which was published 15 years ago: http://eon.sdsu.edu/~seshagir/SDR05.pdf
This filter reportedly results in an order of magnitude lower EVM for the same number of taps, and importantly significantly reduced rejection of adjacent channels, but to my understanding isn't commonly used since the Root-Raised Cosine filter is so baked into our standards.
Is anyone aware of actual use of this alternate pulse shaping filter and does anyone have further experience with it? Are there any other reasons, now 15 years later, that this hasn't been more widely adopted?

Liguid-dsp implements root-Nyquist harris-Moerder (hM-3) filter using Parks-McClellan algorithm as :
void liquid_firdes_hM3(unsigned int _k, unsigned int _m, float _beta, float _dt, float * _h)


Analog Filter Prototype to Direct Form Second Order Cascade I've been using the Faust programming language a lot lately for experimenting with DSP, and I've been digging into their implementation of an Elliptic (Cauer) Lowpass filter. Here's an example of such a filter with order three.

If you're unfamiliar with Faust, hopefully the example there is still fairly clear: they've designed a 3rd order Cauer filter with the properties listed in the comments and have implemented it as a second order direct form filter feeding into a first order direct form filter with the coefficients as listed.

I understand the analog prototype they've designed with [z,p,g] = ncauer(Rp,Rs,3); in Matlab/Octave, and I think I have a good understanding of how to factor the transfer function into cascaded first- and second-order filters. What I don't understand is their use of poly in Octave to find the coefficients there, and why the frequency of this elliptic filter is only governed by the last coefficient in the final first order filter? Intuitively it seems like so many more of those coefficients should depend on the frequency of the filter, especially if we have to consider multiple sample rates. Did they skip some steps here or make some assumptions that I'm missing? Any explanation to my confusion here would be greatly appreciated, thank you!

Update: I find it interesting that in Matlab/Octave I can compute [z,p,g] = ncauer(0.2,60,3); sos=zp2sos(z, p, g), deriving coefficients for second order sections without ever specifying the cutoff frequency. Maybe that's where I'm most confused?

The coefficients in the function lowpass3() are the numerator and denominator coefficients of the two sections (first and second order) of the normalized continuous-time transfer function. The continuous-time transfer function looks like this:

$$H(s)=\frac{b_{21}s^2+b_{11}s+b_{01}}{s^2+a_{11}s+a_{01}}\cdot\frac{1}{s+a_{02}}\tag{1}$$

Note that this normalized filter has  a cut-off frequency of $1$. The two functions tf2s and tf1s transform these two sections to the discrete-time domain using the bilinear transform, whereby the normalized cut-off frequency of the continuous-time filter is mapped to the desired cut-off frequency in the discrete-time domain. So whereas the coefficients in $(1)$ are independent of the cut-off frequency (because the filter is normalized), the coefficients of both final filter sections in the discrete-time domain do depend on the cut-off frequency.

Your example of using ncauer followed by zp2sos will not give you the coefficients of the discrete-time filter with the desired cut-off frequency. Instead, you just the compute the coefficients of the first and second order sections of the normalized continuous-time low pass filter, i.e., the coefficients hard-coded in lowpass3(). What you miss in that case is the transformation step from continuous-time to discrete-time, and the mapping of the normalized cut-off frequency to the desired cut-off frequency.


similar tools for MATLAB ? I have DSP in my academics , I am here to have a suggestion regarding the better alternative to MATLAB because of Its downloading size.

Thank you.

Octave is a good option. Some people call it as a free version of Matlab.
I found Wolfram Mathematica to be of great use. It seems to be much more powerful than Matlab, not to mention, simpler to use.


Why do I hear a "blip" at the end of an audio clip of a sinusoid? I generated a one second audio clip of a sinusoid at 440 Hz (16-bit, 44.1 kHz), ensuring that the last sample is zero. However, I hear a faint "blip" sound at the end of the clip, and have noticed this with similar clips that I have found online. Do you know why this happens, and how I can get rid of the "blip"?

Thanks for your help!

There's still the first derivative discontinuity. That gives a blip with a -12 dB/octave decaying spectrum. It's less audible (because it has less high frequencies) than the -6 dB/octave decaying spectrum of a proper discontinuity, but still audible, especially with a subwoofer...

Do a (long enough a) fade out like Fat32 says and there will be no audible click.


Why do we use Difference RGB to grayscale function? I found this function for RGB to Gray scale converter at this kind:

  RGB to grayscale

  This example converts an image with RGB channels into an image with a
  single grayscale channel.

  The value of each grayscale pixel is calculated as the weighted sum of
  the corresponding red, green and blue pixels as:

  Y = 0.2125 R + 0.7154 G + 0.0721 B

  Copy to clipboard

  These weights are used by CRT phosphors as they better represent human
  perception of red, green and blue than equal weights. 1

Also at Matlab site said about :

  Algorithms

  rgb2gray converts RGB values to grayscale values by forming a weighted
  sum of the R, G, and B components:

  0.2989 * R + 0.5870 * G + 0.1140 * B

  These are the same weights used by the rgb2ntsc function to compute
  the Y component.

  The coefficients used to calculate grayscale values in rgb2gray are
  identical to those used to calculate luminance (E'y) in Rec.ITU-R
  BT.601-7 after rounding to 3 decimal places.

  Rec.ITU-R BT.601-7 calculates E'y using the following formula:

  0.299 * R + 0.587 * G + 0.114 * B

So, i like to know why difference value and how many standard exist in this function type?

Thanks.

There are different color spaces. Each of them has different coefficients for calculating gamma-corrected luma from gamma-corrected RGB values. You have listed coefficients for calculating luma according to ITU-R BT.709, ITU-R BT.601 and NTSC color spaces. There are many others as well, like BT.2020 used in consumer TVs for Wide Color Gamut.


Can we generate signals in reality that have non symmetric magnitude spectrum in the Fourier domain? I am reading communication systems, and have a doubt:
Is it possible to generate signals in reality that have a non-symmetric magnitude spectrum in the Fourier domain?
For example, if I have a signal
$$u(f)$$ or something like this as Fourier domain representation
I know that their Inverse Fourier transform is going to have imaginary element $j$ in them, but does this $j$ just symbolize 90-degree phaseshift, can I generate signals in reality that have asymmetric magnitudes of their Fourier transforms.
I have read that $j/-j$ just represent a 90-degree lag or lead, but do such signals with asymmetric magnitude spectrum in the Fourier domain exist in reality?
This will also clear my doubt that whether there are such signals in reality that have just positive/negative frequency components in their Fourier domain representation.
As the answers suggest having $j$ just means phase shift, but if some practically permissible signal has representation $jm(t)$, where $m(t)$ is real, with no imaginary component, then with respect to what I have to define the phase.
Or such a condition will never arise where the inverse Fourier transform of a signal is purely an imaginary quantity.

As @Bob said, it depends of your definition of "in reality". I am going to ignore images on this post and focus on signals as a function of time.

One-dimensional signals always have a symmetrical magnitude spectrum. Examples are the voltage or current on a wire, or a sequence of real numbers stored in digital memory.

Two-dimensional signals may have a non-symmetrical magnitude spectrum. Examples are complex signals, whether defined mathematically or stored in a digital memory; and the voltage or current on two wires, one of which is labeled $\sqrt{-1}$ and with the agreement of doing complex arithmetic (see $9\rm V$ Battery with $45^\circ$ phase).

In my opinon, all these signals exist in reality, and then the answer is that yes, some signals in reality may indeed have a non-symmetrical magnitude spectrum.


Fourier Transform of a signal using direct integration and properties Am trying to compute the Fourier Transform of a function using the properties of the Fourier Transform once and checking my answer using direct integration. My problem is that am not getting the same for both ways.

Sorry for not being able to type the equations as I am not familiar on how to do that, It would be great if somebody can tell me how. Your help is appreciated. Thanks.

Your first solution using the properties of the Fourier transform is correct. Your second solution is wrong, because you forgot to include the unit step function. Your function $g(t)$ should be defined by

$$g(t)=e^{-t}u(t)\tag{1}$$

which gives for $g(2t-1)$

$$g(2t-1)=e^{-(2t-1)}u(2t-1)=e^{-(2t-1)}u\left(t-\frac12\right)\tag{2}$$

Consequently, the Fourier integral is given by

$$\mathcal{F}\{g(2t-1)\}=\int_{\frac12}^{\infty}e^{-(2t-1)}e^{-j\omega t}dt\tag{3}$$

I'm sure you can prove now that the result of $(3)$ equals the result you got from using the properties of the Fourier transform.


Optimized 2D wavelet transform using FFT I'm currenty aiming to optimize my fast wavelet transform (FWT) algorithm for 2D signals (images). It works as follows:

one iteration of 1D FWT does convolution of 1D input data with a selected 1D filter (lengths from 2 to approx. 60) and downsamples the result
algorithm for 2D transform does 1D FWT across all rows and then all columns of the input image
it iterates if more levels are desired

The transform is a part of interactive application that demonstrates wavelets and their use. It works fairly fast and usually responds in real time to user's interactions. But if the filter is very long some performance issues occur. I've read that using Fast Fourier Transorm (FFT) instead of convolution is effective for long enough filters.

I've already implemented the 1D FFT, but the question is how to use it for maximum efficiency? Should I transform the input data before single 1D FWT, then perform convolution (which corresponds to multiplication in frequency domain), and then transform the data back using inverse FFT? Also, how is the multiplication done exactly? For example, the input data of length 256 and filter of length 4 are both transformed using FFT and then only the first 4 values of input data are multiplied before transforming the data back? I'm struggling a little bit on the details and would very much appreciate any insight into this.

EDIT:
I've figured out that in my case I'm after circular convolution, therefore the filter should be zero padded so that its length is the same as length of the input data. But my question about efficiency still holds. How should I use FFT for FWT computation in order to be beneficial?

I'm still not convinced you are spending your time optimally. Probably you would get better gains doing other optimizations. Already 10-15 taps are considered very long wavelet filters for most applications that I have seen. But nevermind that I will try and answer your question.

You can do two different approaches

You could do it one level at a time.
Everything at once.

There are drawbacks and gains from either of them. I will start with the first one.

If you go for one level at a time:

You zero pad filters to be equal lenght as the rows and columns. Also make sure you pad so the filter gets in the middle.
FFT the padded filters and corresponding image dimension. If width $\neq$ height you would need to do 2 FFTs per filter one of each dimension length.
Pointwise multiplication of the FFTs of filter and image.
Inverse FFT of the product image.
Downsample product image in the dimension you transformed.
Repeat for the next dimension.
Repeat for all levels.

Notice that this approach will involve a set of iFFTs and FFTs at each level. That would be double-bad.

The everything-at-once approach:

Same as above but instead of iFFT and decimating each level to just FFT back again on the decimated signal you just multiply several times with the filter and do iFFT and decimate 2,4,8 times et-cetera. Then you only need to do one set of FFTs (at the start) but still a full set of iFFTs for each level.

Also the multiplications while in the Fourier domain will be more computationally expensive than ordinary multiplications as the fourier components will be $\in \mathbb{C}$ but filter and image coefficients $\in \mathbb{R}$ ( I guess ) and you would also likely need higher numerical precision to get the same error.


How do I time-stretch a signal from the frequency domain? I have $FFT(signal)$, and want to find $operation$ such that $$IFFT(operation(FFT(signal),2)) = timestretch(signal, 2)$$
Where $timestretch$ stretches the signal in the time domain - maintaining the frequencies, unlike basic subsampling stretching which leads pitch change.
One option I know I could do is convert to a 2D spectrogram and stretch this "image" and then reverse the spectrogram, but this seems both inaccurate and expensive.

Are you talking about taking an audio recording where someone was speaking too fast, and you want to slow down the speech without changing the pitch? If so, then your original formulation of the problem is incorrect; you can’t do what you want by taking one long FFT of the entire recording and then processIng it.  The usual way this is solved is by using the phase vocoder, where you break the signal into overlapping blocks and transform each block.  This is a deep subject with lots of tricks to avoid artifacts.


comparison of speech signal processing and financial data I have read that in speech signal processing analysis when voice is segmented in brief temporary segments the series segments transitions from being non stationary to stationary. My question is if this is also the case for financial data when this data can be segmented in smaller and smaller time spand windows (stock market).

Financial market trading data is thought to be fractal in dimension and thus not sufficiently band-limited before sampling.  So smaller windows might not have any effect on improving pseudo-stationarity.

Whereas the physical processes behind speech production are better understood and more physically limited (the rate of human muscle contractions controlling lung, larynx, and head cavity resonances, etc.), as well as the sound capture being low-pass filtered below typical human hearing thresholds before sampling above the Nyquist rate.


How do I play and record audio at the same time? So this question is more related to Matlab itself.

Some background info: I'm trying out the Exponential Sine Sweep method to obtain the Room Impulse Response with my laptop speakers and microphone for testing.

I can generate the Sine sweep perfectly but I have trouble playing and recording it without a lot of lag. I'm now trying to use "audioPlayerRecorder" from the Audio System Toolbox but when I do:

f1 = 80;       % Start frequency (Hz)
f2 = 22000;    % Stop frequency (Hz)
T = 3;         % Sweep duration in seconds
fs = 44100;    % Sampling frequency

[sweep, inverse] = ESS_generator(f1, f2, T, fs);

playRec = audioPlayerRecorder(fs);
recording = playRec(sweep);

I get an error on the last line saying:
Error using coder.internal.assert (line 33)
No full-duplex audio device detected

But I know that my laptop is able to use the microphone and speakers simultaneously so I'm thinking maybe the error gets thrown because of a missing configuration?

I'm running the code on Linux with the ALSA audio drivers. I'm using external speakers through the AUX audio port, and the default laptop microphone.

I suspect that your problem is to do with having your speakers and microphone on different ALSA audio devices or subdevices.

You need to find the available sound cards for recording and playback on your system. The following commands will tell you this information :

aplay -l
arecord -l

You will need to define a default full duplex sound card using your choice of devices you see in the output to the above commands.

The way to define this full duplex card (made of two different soundcard devices) is to define it in your ~/.asoundrc file. This example shows you how to do that, however I will repeat it here.

The requirement is that you use plughw as your device name because the clocks on the two sound cards will drift and you will need resampling of some form. But you can play with this option.

So for example if you want to take input from a device on the first card (hw:0) and output to the second device on the second card (hw:1,1) then you would write the following in your ~/.asoundrc file :

pcm.!default {
    type asym
    playback.pcm "plughw:1,1"
    capture.pcm  "plughw:0"
}
ctl.!default {
        type hw
        card 1
}

Hopefully now Matlab sees a full duplex card as the default device and allows full duplex operation.


What is the relationship between frequency vector and the peak frequency in FFT? I am fairly new to both Matlab and FFT. I would like someone to explain to me the relationship between the frequency of the wave being processed, the number of point (usually called NFFT), the Nyquist Fs and frequency vector. Basically how do you find the peak frequencies? how does FFT represent the frequencies? Any and all help is welcome. Thank you

Although I think Jason's answer contains all the information you requested, maybe what you're looking for is a little less theoretical.

Let's assume you have a signal of $$x(t)=cos(2\pi t f_x)$$
where $t$ is the time in seconds and its frequency $f_x=20$Hz. Let's assume further you sample this signal at $f_s=1000$Hz over $T=3$ seconds, i.e. you have $K=3000$ samples $x[k]$, where $k$ is the integer sample index of the signal $x(t)$.

The sampling theorem basically says that your sampling frequency $f_s$ must be at least twice the maximum signal frequency (which is the Nyquist frequency $f_s/2$) in order to reproduce it without aliasing. For our example: $f_s=1000>2f_x=40$.

You choose the "number of points" $N$ for the FFT (you called it NFFT) yourself, depending on how many different frequency samples you need to distinguish. The maximum frequency resolution you get by using all available samples you have for calculating the FFT, which - depending on the $f_s$ and $T$ of your signal - might not be feasible.

For our example, using all 3000 samples $x[k]$ is not a problem on a PC running Matlab, so let's assume we do that get $N=3000$ complex frequency samples $X[n]$ where $n$ is the integer frequency index.

Now to your question where in $X[n]$ to find a specific frequency $f$ you're interested in if $N=K$: $$n = f T$$

Therefore, the Nyquist frequency $f_s/2$ bin is at $n_s = T f_s/2 = 1500$, right in the middle of the FFT result vector. If you plot the amplitude spectrum $abs(X)/N$, you'll notice that there are peaks at $n_x = T f_x = 60$ and $n_x' = 2940 = N - 60$, each of magnitude $1/2$. This is because the spectrum is "mirrored" around $f_s$ and the signal amplitude of 1 is divided between both peaks $n_x$ and $n_x'$.

For Matlab indices, I think you'll have to add 1 to the $n$ mentioned above because Matlab's vector indices start at 1 instead of 0.

Note that if you have frequencies $f$ which are not a multiple of your sampling frequency $f_s$, the amplitude will spread to neighbouring bins and you'll have to round $n$ to the nearest integer.

See also Matlab's FFT documentation, which has a nice example.


Calculating phase shift between two sines of equal, unknown, amplitude and frequency I've got a source signal $s_1(t)$ and a system affecting that, yielding a phase-shifted version $s_2(t)$;

\begin{align}
\color{blue}{s_1(t)} &= A\sin(t+\phi_1)\\
\color{red}{s_2(t)} &= A\sin(t+\phi_2)\\
\end{align}

The amplitude $A$ is unknown. How many samples (sampling happens periodically) would I need to verify that $\phi_1-\phi_2 = \frac32\pi$?

As a practical example, given these 2 points on each curve: on the purple curve, the original signal:

point A (231,3.533151)
point B: (230,3.614015)

and on the the filtered signal: the red curve:

point C (231,1.398873)
point D (230,1.174265)

How would you confirm the phase shift of $\frac32\pi$, given an equal -unknown- amplitude and equal -unknown- frequency of the 2 signals? Is a numerical solution possible for online, causal analysis? Or do you need 3 points on each curve? Thank you.

Two samples is sufficient and we can determine if the phase difference is indeed $\frac{3\pi}{2}$ as follows:

Start with the hypothesis that the phase difference $\phi_1-\phi_2$ is $\frac{3\pi}{2}$, which is equivalent to $-\frac{\pi}{2}$

If and only if the phase is in such quadrature, then $s_1(t)$ and $s_2(t)$ will be the real and imaginary components of a single vector on the unit circle with radius A which can be given as $-jAe^{j(t+\phi_1)} = As_1(t) + jAs_2(t)$. This  starts at the angle $\phi_1-\pi/2$ when t = 0 and rotates counter-clock-wise with increasing t, which matches the plot shown with the purple line representing the real axis and the red line representing the imaginary.

Starting from a position on the unit circle at time = $t_1$ with position given by the real and imaginary components $s_1(t_1)$ and $s_2(t_1)$, we can move forward in time to $t_2$ over a determined angle which should predict $s_1(t_2)$ and $s_2(t_2)$. If this prediction matches the given result, then this will confirm the hypothesis:

The starting angle using the first sample as given by B and D is given by $tan^{-1}(s_1(t_1)/s_2(t_1))$
$= tan^{-1}(3.614015/1.174265) = 1.256637 $ radians.

The change is angle to the second sample is simply the difference in time since the frequency unit is 1 radian/sec given by $sin(t)$:

$\Delta\phi = t_2-t_1 = 231-230 = 1$ radian    (NOTE! This alone clearly does not match the plot visually, so the numbers given cannot be actually derived from the plot! There must be a frequency factor multiplied by t in the actual formulas for this to match the plot).

Therefore if the phase difference was $-\frac{\pi}{2}$ then the second angle would be at $1.256637 + 1$ radian or $2.256637$ radians and the ratio of $(s_1(t_2)/s_2(t_2))$ would be $\tan(2.256637)= -1.22195$

This does not match the ratio for the second sample given by A and C and therefore the phase difference cannot be $\frac{3\pi}{2}$ for the formulas as given.


What is the difference between transport delay and dead time? I am confused on what the difference between transport delay and dead time is. They seem to have the same function, but I'm not sure. Can someone explain the difference if there is one?

According to this document (p. 4) and this one (p. 13), transport delay and dead time denote the same thing: latency.


How to detect continuous noise in audio call? We want to detect bad audio quality in audio call.

Should we extract this information from amplitude or frequency spectrum ?
Is there any technique/transform/filter available for this ?

Example:

I am attaching two signal graph(Clean Signal vs Noisy One) of Amplitude vs time.

Clean Signal Graph:

Noisy Signal Graph:

As shown in above graph, there is a continuous noise in background in Noisy signal.

How can it be differentiated between this two signal ?

My recommendation would be to do the processing in the frequency domain as methods are available that can directly be used to approach your problem. In many cases, speech "quality" is related to the signal-to-noise ratio (SNR) of the signal. At least, I assume this is the case for your application - i.e., telephone communication. In general, some measure for the "quality" can be derived from the SNR then with low SNR indicating "bad quality" and, conversely, high SNR indicating "good quality".

One approach to solve your problem would hence be to try and estimate SNR of the recording. In order to estimate the SNR, well-known noise power spectral density (PSD) estimation algorithms such as, e.g., the "minimum statistics" algorithm [1] or a minimum mean-square error based algorithm [2] can be used. An estimate for the SNR can then be computed from the noise PSD as

$$
\hat{\text{SNR}} = \frac{P_{xx}}{\hat{P}_{nn}}-1,
$$

with $P_{xx}$ the PSD of the noisy input signal and $\hat{P}_{nn}$ an estimate of the noise PSD. A single SNR estimate can be computed for the complete signal, but as the SNR changes strongly over time, it might be useful to compute an SNR estimate for short blocks of the signal - typical block lengths are in the order of $20\,$-$30\,\text{ms}$.

Instead of first estimating the noise PSD, the SNR can be estimated directly, e.g., with the method proposed in [3].

Alternatively, some non-intrusive (i.e., no clean reference signal is required) measure for either the speech quality or the speech intelligiblity can be used. You could try the speech-to-reverberation modulation energy ratio measure [4], of which implementations in Python and MATLAB are directly available.

References

[1] Rainer Martin, “Noise Power Spectral Density Estimation Based on Optimal Smoothing and Minimum Statistics,” IEEE Transactions on Speech and Audio Processing, vol. 9, no. 5, pp. 504–512, Jul. 2001. PDF MATLAB implementation

[2] T. Gerkmann and R. C. Hendriks, “Unbiased MMSE-Based Noise Power Estimation With Low Complexity and Low Tracking Delay,” Audio, Speech, and Language Processing, IEEE Transactions on, vol. 20, no. 4, pp. 1383–1393, May 2012. PDF MATLAB implementation

[3] E. Nemer, R. Goubran, and S. Mahmoud, “SNR Estimation of Speech Signals using Subbands and Fourth-Order Statistics,” IEEE Signal Processing Letters, vol. 6, no. 7, pp. 171–174, Jul. 1999. IEEE Xplore

[4] T. H. Falk, C. Zheng, and W. Y. Chan, “A Non-Intrusive Quality and Intelligibility Measure of Reverberant and Dereverberated Speech,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 18, no. 7, pp. 1766–1774, Sep. 2010. IEEE Xplore Python implementation MATLAB implementation


Compressive Sensing vs. Sparse Coding There apparently are different terminologies used to refer to the same field called "compressive sensing" such as (see this wiki page): compressed sensing, compressive sampling, or sparse sampling. I wonder about "sparse sensing" though!

Nonetheless, and after some internet search, what people refer to as "sparse coding" seems to not refer to the "compressive sensing" field as the other terminologies I cited above.

Is there really a difference between compressive sensing and sparse coding?

What about dictionary learning?

As you correctly noted compressed sensing, compressive sampling, sparse sampling all mean the same thing. Some authors also call it sparse sensing.
The idea behind compressed sensing is that a sparse signal can be recovered from very few linear measurements. In symbols, if $\mathbf x$ is $N\times 1$ sparse$^\ddagger$ vector, and $\mathbf A$ is an $M\times N$ matrix with $M\ll N$, and we measure $\mathbf y=\mathbf{Ax}$, then compressed sensing theory tells us$^\dagger$ that we can exactly recover $\mathbf x$ from $\mathbf y$. This is remarkable because it says that we can recover the original signal from fewer measurements.

Dictionary learning on the other hand deals with an entirely different problem of representing a bunch of data vectors in a parsimonious way. Given a set of data vectors $\{\mathbf x_1, \mathbf x_2,\ldots, \mathbf x_K\}$, we would like to find another set of vectors $\{\mathbf v_1, \mathbf v_2,\ldots, \mathbf v_L\}$ (called "atoms") such that that each data vector $\mathbf x_i$ can be represented as a linear combination of these $\mathbf v_j$'s. The set of atoms is called a dictionary. The goal here is to learn a dictionary that is much smaller than the number of data vectors$^*$ i.e. $L < K$.

Given a set of atoms in a dictionary and a vector $\mathbf y$, the goal of sparse coding is to represent $\mathbf y$ as a linear combination of as few atoms as possible.

Finally, sparse dictionary learning is a combination of dictionary learning and sparse coding. The goal here is two-fold: finding a parsimonious representation of the set of data vectors and ensuring that each data vector can be written as a linear combination of as few of the atoms as possible.

Compressed Sensing v/s Sparse Coding
Both of these techniques deal with finding a sparse representation but there are subtle differences.

Compressed sensing deals specifically with the problem of solving an underdetermined system of linear equations i.e. fewer data points than the original signal. From an unknown sparse signal $\mathbf x$ and sensing matrix $\mathbf A$, we observe the data vector $\mathbf y = \mathbf{Ax}$. $\mathbf A$ has fewer rows than columns. Compressed sensing theory deals with the following kinds of questions questions:

Under what conditions is the under-determined set of linear equations solvable and how do we solve it in a noise-robust, computationally tractable manner?
How do we design sensing matrices $\mathbf A$ for various applications?

In contrast, sparse coding does not deal with the question of designing $\mathbf A$. Moreover you aren't interested in solving under-determined system of equations --- $\mathbf A$ is allowed to have more rows than columns.$^\%$

References:

Compressive Sensing [Lecture Notes]

Dictionary Learning

Online dictionary learning for sparse coding

Footnotes:

$^\ddagger$ Sparse means the vector has very few non-zero elements.

$^\dagger$ $\mathbf A$ and $M$ need to satisfy some technical conditions.

$^*$ Unlike standard transform methods such as Fourier transform, dictionary learning is data-adaptive. When taking a Fourier transform, the basis vectors $\mathbf v_j$'s are fixed ahead of time (complex exponentials). In dictionary learning, they are learned from data.

$^\%$ This is called an over-complete dictionary.


Measure how Straight/Smooth the Borders are Rendered in an Image I have two images:

I want to measure how straight/smooth the text borders are rendered.

First image is rendered perfectly straight, so it deserves a quality measure 1. On the other hand, the second image is rendered with a lot of varied curves (rough in a way) that is why it deserves a quality measure less than 1. How will I measure it using image processing or any Python function or any function written in other languages?

Clarification :

There are font styles that are rendered originally with straight strokes but there are also font styles that are rendered smoothly just like the cursive font styles. What I'm really after is to differentiate the text border surface roughness of the characters by giving it a quality measure.

I want to measure how straight/smooth the text borders are rendered in an image. Inversely, it can also be said that I want to measure how rough the text borders are rendered in an image.

I'd try a very "tinkery" approach here:

Erode the image, so that the black area is shrunk by a fixed radius of pixels from its border (say, 5px).
Dilate the resulting image by the same amount
measure the amount of difference between original and processed image.

The idea is that something that is a locally convex border doesn't suffer through erosion (it's only shrunk) significantly, and that this erosion can be reverted by dilatation.


Speech recognition and comparison algorithms I am doing several research for a prototype which I will be building. My prototype will consist of single word speech recognition and comparison with a preset recorded words. The prototyped system will then ask the user to say a the word and give a percentage match.

I am hoping to build this prototype with JAVA, however I am not at that stage yet, I am still researching all the available algorithms that exist in order to implement the above. I am seeking some guidance from you experts out there if I am on the right track in looking at the following algorithms to be able to implement the above:

Acoustic speech recognition
Hidden Markov Model algorithm
Discrete Fourier Transform
Fast Fourier Transform
The Cross-correlation Algorithm
The Auto-correlation Algorithm
The FIR Wiener Filter
Dynamic Time Warping algorithms
Mel Frequency Cepstral Coefficient

I am not sure if the above is wrong for my project and look into audio fingerprinting.

I am sorry I cannot provide any code for the time being I am still in the theory research part, If this question is not suitable to this forum please let me know so I can delete it and waste no one's time as it is not my intention.

P.S. Keep in mind that since this is only a prototype, moving forward I will research and eventually use any already existing algorithms (API's/Libraries in JAVA) so that is mostly why I am asking if I am on the right track. Many thanks for your help guys !

What you want to do sounds like standard speech recognition. I recommend you check out the Open Source JAVA library from CMU, Sphinx4. Run the Transcriber or Dialog demos and you can see it in action doing some basic recognition of digits & simple english commands from a small grammar. Note that many of your algorithms listed above, fit into different parts of a speech recognition system (frontend processing, search manager, acoustic scorer, etc)

Check out this paper for a good overview (and to see how all your algorithms above relate)


Possibilties of chebyshev polynomial waveshaping When talking about harmonic distortion, or more specifically waveshaping, we say the order of distortion, can be solved from the equation:

$x^n$

Where n is the order of distortion. To my knowledge the complete result of any waveshaping function can be had by summing all of the different orders of harmonic distortion, multiplied with an appropriate index $a_n$.

$$\sum_{n=1}^{\infty} a_nx^n = WaveshaperFunction$$

Be aware that by feeding the equation with $x=sin(\omega t)$, the highest frequency any $x^n$ can produce is $sin(n*\omega t) $. However, higher orders can still produce multiple different frequencies for a sinusoidal input.

Chebyshev polynomials addresses the issue, removing any other frequencies limiting the harmonic distortion to only $sin(n*\omega t)$, producing an orthogonal result. The polynomials are defined as:

$T_0(x) = 1$

$T_1(x) = 2x-1$

$T_{n+1}(x)=2xT_n-T_{n-1}$

Now, since the result of putting any $sin(\omega t)$ through the chebyshev waveshaper is a sine with a frequency multiple of N, it is easy to say that any harmonic signal with the original signal's phase (or inverted phase) can be produced by summing over all the chebyshev polynomials. That is to say, summing over all the polynomials, any waveshaping function can be generated.

$$\sum_{n=1}^{\infty} b_nT_n(x) = WaveshaperFunction$$

Where $b_n$ represents the amplitude of nth generated harmonic, given a sine input. $T_n$ represents the cheyshev polynomial producing nth harmonic given sine input. I think there could be many advantages for thinking about time-invariant harmonic distortion or intermodulation distortion through this framework rather than through the before-mentioned orders.

As for the question, did I make a mistake somewhere? Particularly I am not so certain that any function can be generated through summing the chebyshev polynomials. Chebyshev polynomials seem to not care about the amplitude of the input sine wave (apart from some DC offset differences), yet I can easily generate a waveshaper function that dose so (any function that is linear near 0 and then changes). Re-framed, can the polynomials generate any possible function?

For $-1 <= x <= 1$, let's compare Chebyshev polynomials of the first kind, $T_n(x)$, and the basis functions of the Fourier cosine series, $F_n(x)$:

 $F_n(x)=\cos(n \pi x)$

 $T_n(x)=\cos(n\ \text{acos}\ x)$

Writing $T_n(x_T) = F_n(x_F)$ and solving for $x_F$ gives $x_F = (\text{acos}\ x_T) / pi$, revealing that the Chebyshev polynomial series is simply an argument-warped version of the Fourier cosine series that maps the range $0 <= x_F <= 1$ to $-1 <= x_T <= 1$:

This exclusion of $x_F < 0$ makes irrelevant the constraint that the Fourier cosine series is only applicable to even functions, functions symmetric around $x_F = 0$, meaning that $f(x_F) = f(-x_F)$. The mapping function is well-behaved in the range of interest, so if the Fourier cosine series can represent any even function, then the Chebyshev polynomial series can represent any function.

Not knowing the amplitude of the input sinusoid, it is not possible to control the proportional amplitudes of the generated harmonics by the choice of the waveshaper. For example if the waveshaper is $f(x) = x^2 + x$:

$$f(\cos x) = \frac{1}{2} + \cos x + \frac{1}{2} \cos 2x\\
f(\frac{1}{2} \cos x) = \frac{1}{8} + \frac{1}{2} \cos x + \frac{1}{8} \cos 2x$$

The first one has the amplitude of the second harmonic at -6 dB compared to the fundamental, and the second one at -12 dB. Chebyshev polynomials are no miracle cure. The above polynomial $x^2 + x$ was a weighted sum of the first three Chebyshev polynomials. Any polynomial can be written as a weighted sum of Chebyshev polynomials, and vice versa.

The problem persists for pure Chebyshev polynomials, for example $f(x) = 4 x^3 - 3x$:

$$f(\cos x) = \cos 3x\\
f(\frac{1}{2} \cos x) = \frac{1}{8} \cos 3x - \frac{9}{8} \cos x$$

A Chebyshev polynomial series does have the advantage of better numerical stability than a simple power series, and if the amplitude of the input sinusoid is controlled somehow, it can be used to control the level of harmonics in the way you describe.


Image reconstruction using Irradiation I am a bit confused as to how to go about this.

I have a surface [x,y,z] given by:

[x, y] = meshgrid(-10:10,-10:10);
z = x .* exp(-x.^2 - y.^2);
surf(x, y, z);

I want to use $I=\rho n^TS$ where $\rho=0.1$, $S$ is the sources vector and $n$ is the normal to the surface, to reconstruct the image. Say I want to do this for 10 arbitrary $S$ vectors and reconstruct 10 images. how do I do this?

I tried extracting $n$ by [nx, ny, nz] =surfnorm(x,y,z); but I am confused here as to how to represent the $S$ vector. is it supposed to $2\times 21$? for one vector? Is there an easier way I am missing?

Each point in surface has a normal represented by 3x1 vector so transpose of it is 1x3 . Considering radiance from the source is represented by normal to the source it will be represented by a 3x1 vector.


Intuitive meaning of complex time domain signal representation (quadrature sampling) I saw What is the meaning of imaginary values in the time domain?.

My question is about I Q samples as recorded from an SDR receiver (e.g. RTL-SDR).

My understanding is that we can call these I Q samples a complex time domain signal and these can have both non-zero real and imaginary components (the question I linked above is talking about zero imaginary part in the time domain).

Is there an intuitive explanation of what is the meaning of this type of complex time domain signal?

I can think of taking a simple case of an input sine wave of a frequency below the sampling frequency which would make the IQ signal's phasor rotate with the frequency equal to the difference between the input sine wave's frequency and the sampling frequency thus allowing to unambiguously reconstruct frequencies up to (but not including) the sampling frequency from the IQ samples.

But I'm struggling to extend this special case to an input composed of several sinusoids of different frequencies.

One common use of the IQ data from an SDR is to feed them to an FFT.  The upper and lower halves of an FFT result indicated different relationship between the (real) cosine and (imaginary) sine components of the complex FFT result. The direction of rotation of each complex component tells the rest of the signal chain whether the complex component is an upper or lower sideband of the assumed center tuned frequency of the IQ SDR data stream.

The quadrature component in the IQ input is what allows the two halves of an FFT result to not just be redundant complex conjugates of each other (which they would be given a strictly real input, no Q or zero Q).  Thus allowing separating the upper and lower sidebands as independent spectrum.

Also see: https://electronics.stackexchange.com/questions/39796/can-somebody-explain-what-iq-quadrature-means-in-terms-of-sdr/99617#99617


Why would super resolution using deep learning beat the old school techniques? I assume a standard/ old school super resolution techniques use Polynomial interpolation  and Upsampling approaches. Why would deep neural networks achieve better results?

The objective is to derive a higher resolution image from a low resolution image. When performing this upsampling, inevitably, you have to fill in information that you do not have.

The differences between the two approaches are in how you fill that information in.

In the standard techniques an underlying model is assumed. So, for example, in nearest neighbour interpolation we decide to fill in the unknown pixels with the value of the nearest pixel. An improvement over that is linear interpolation where we derive a weighted sum of the pixel values around the unknown pixel to set its value. An improvement over that are spline and cubic interpolation which are trying to derive the value of the unknown pixel from a bigger patch of known pixels around it.

In all of these cases, we are assuming a model for filling in the pixel values which is based on physics (e.g. due to the Point Spread Function) or assumptions.

But, how about deriving the "model" that constraints the values of the pixels from the images themselves and then use that model and the limited known data to guess the unknown pixels? Instead of deciding on an interpolant arbitrarily, we derive it.

This is an optimisation problem where we are trying to discover the function that connects the low and high resolution image versions. To do this, you obviously need both the high and low resolution images. But what you can do, is take (known) high resolution images, degrade them (reduce their resolution) and workout the relationship between the high resolution and low resolution.

The surprising result here is that images depicting different content end up having very similar interpolants. Of course, content does play a role in the sense that if you derive models from a set of images with a specific bandwidth, you will discover models that support that specific bandwidth. You cannot train the model on smoothly varying images with lots of low frequency content and expect it to reconstruct high contrast high details parts of unknown images.

Therefore, the technique is more suitable to classes of images and so far it has been applied to medical imaging (and of course Magnetic Resonance Imaging where there is a long history of reconstruction from limited data) and radio astronomy

For more information you might want to see this and this.

Hope this helps.


Why idft(dft(a) * dft(b)) not equal to convolve(a, b)? I'm a little confused... I always thought the DFT of a convolution was equal to a product of DFTs, but when I tried this in Python:

from scipy import *

a = [1+0j, 2+0j]
b = [4+0j, 5+0j]

print list(ifft(fft(a) * fft(b)))
print list(convolve(a, b))

I got back:

[(14+0j), (13+0j)]
[(4+0j), (13+0j), (10+0j)]

Why are they not the same thing?

ifft(fft(a) * fft(b)) performs a cyclic convolution, convolve apparently zero-pads the inputs. If you pad both arrays with zeros, the result should be the same:

a = [0,0,0,1+0j, 2+0j,0,0,0]
b = [0,0,0,4+0j, 5+0j,0,0,0]

print list(ifft(fft(a) * fft(b)))
print list(convolve(a, b))


How to find lattice coefficients of all pole IIR filter with filter coefficient $k_m$ is 1? Given the following  system function
$$ H(z) =\frac{1}{1+\frac{2}{3}z^{-1}+\frac{5}{8}z^{-2}+\frac{2}{3}z^{-3}+z^{-4}} \tag{1} $$
We have to draw lattice structure for above system function

The denominator polynomial of the transfer function at lattice stage number $m=4$ is assumed as $A_m$  and now $A_{m-1}$ can be found with

$$\begin{align} A_{m-1}(z)&=\frac{A_{m}(z)-k_mB_{m}(z)}{1-k_m^{2}}\\
\text{with}\quad B_m(z)&=z^{-m}A_m(1/z)\\
\text{and}\quad A_0(z)&=B_0(z)=1\end{align}\tag{2}$$

where $k_m$ is the $m^{th}$ reflection coefficient.
But since the $k_m$ value is 1 equation (2) becomes zero in denominator and hence cannot be solved.

I couldn't find the solution to this particular problem in some of my book or other site.
How can I get the lattice coefficient for above system function?


Accurate description of the spectrum of sinusoid waveform added with WGN? The transmitting signal is
$$
s(t) = \sin(20\pi t)+\sin(40\pi t)+\sin(60\pi t)
$$

The received signal after AWGN channel is
$$
r(t) = s(t)+n(t)
$$
where the variance of noise $n(t)$ is $\sigma^2$.

The theoretical single  spectrum of the $r(t)$ should be

$$
R(f) = \pi\delta(f-10)+ \pi\delta(f-20)+ \pi\delta(f-30)+\sigma^2
$$

where $\delta(.)$ is Delta function. However, the spectrum obtained through MATLAB seems not like the theroetical one.

When the SNR is setted as lower than -25dB, the spikes will be submerged by noise. And how to calculate the threshold SNR where sinusoid part is totally covered by noise?

The relative code is here

fs_init = 10; % basis frequency
om=fs_init*2*pi;
T = 2*pi/om;
tall = 10;
h=T*0.01;
t=0:h:tall-h;
fs_init= 10;
ut = zeros(1,length(t));

%%% Added constants
NSyboml = 3;
a=1;

for i = 1:NSyboml
    ut = ut+a*sin(2*pi*i*fs_init*t);
end

SNR = 0;
rt = awgn(ut,SNR,'measured');   % add noise

% FFT
y = rt;
L = length(ut);    % Length of signal
Ts_v = t./(0:L-1); % sample time
Fs = ceil(1/Ts_v(2));    % sample frequency
deltaf = 0.1; % resolution
NFFT = floor(Fs/deltaf);
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,floor(NFFT/2+1));

% Plot single-sided amplitude spectrum.
 figure;set(gcf,'Color','w','Position',[0 0 500 350]);
plot(f,20*log10(abs(Y(1:NFFT/2+1))));
xlim([0,6*fs_init])
xlabel('Frequency (Hz)');ylabel('|R(f)|');
grid on;
ylim([-120,0])
title(['SNR = ',num2str(SNR),'dB'])
%%%enter preformatted text here

Notice that in your first image, there are peaks where they would be expected to be, so that's fine. I guess you think that "the spectrum obtained through MATLAB seems not like the theroetical one" because of the noisy spectrum you get, instead of a straight line at $\sigma^2$ on the vertical axis.

The random signal added to $s(t)$ is not random actually, it's just a realization you get when calling awgn(). So the FFT will return just a noisy spectrum for that signal. If you generate lots of realizations of $r(t)$ and plot the average of all their FFTs, you'll obtain what you were expecting for.

For your second question on how to calculate the threshold SNR, some criteria is needed. If you average lots of FFTs the variance of the noisy part of the spectrum will drastically decrease, so the SNR needed to see the peaks will be much smaller. If you use just one realization, then the SNR will be higher. There is not a fixed number for the needed SNR, it depends on what you refer to "sinusoid covered by noise". If the spectrum observed was just like the theoretical one (i.e. you averaged a lot of FFTs), a sinusoid would be "completely covered by noise" when its amplitude is such that the corresponding peak in the frequency domain is below $\sigma^2$.


Hilbert Transform of an FM Signal The Hilbert Transform of an Amplitude Modulated signal returns the envelope of the signal. What does the Hilbert transform of a Frequency Modulated signal return? How can I use the Hilbert Transform to get the sidebands of a Frequency Modulated signal?

From this great answer . .

  If we have a real-valued sequence x(n)x(n), and the Hilbert transform
  of x(n)x(n) is xH(n)xH(n), the analytic signal, a(n)a(n), associated
  with x(n)x(n) is:

  a(n)=x(n)+jxH(n)

  Computing the magnitude of a(n)a(n) will give you the instantaneous
  envelope of the original x(n)x(n) sequence. This works great for AM
  demodulation. If you compute the instantaneous phase of a(n)a(n) to
  get ϕ(n)ϕ(n), and compute the time derivative of ϕ(n)ϕ(n), you have FM
  demodulation.

A practical example where we generate an FM signal in Matlab and analyse it looks as follows:

fs = 44.1e3; %// sample rate
dt = 1/fs;
fc = 200; %// carrier
df = 50; %// modulation excursion (Hz)
fm = 4; %// modulation rate (Hz)

tAx = dt:dt:1; %// time axis in seconds
u = sin(2*pi*fc*tAx + (df/fm)*cos(2*pi*fm*tAx));
w = unwrap(angle(hilbert(u)));
v = diff(w)/(2*pi*dt); %// instantaneous phase differential normalized to Hz

Plotting u and v then gives the following plot, where the sidebands occur at the extremes of the signal . .


Expectation of product of hermitian matrix Let $\textbf{h} \in \mathbb{C}^{N\times 1}$, $\textbf{a} \in \mathbb{C}^{N\times 1}$, $\textbf{b} \in \mathbb{C}^{N\times 1}$ and $\textbf{c} \in \mathbb{C}^{N\times 1}$. The variable $h_i$ is defined as $h_i=\left(a_i+b_i\right)c_i$

$\textbf{a}=[1~2~3~4]^T$, $\textbf{b}=[5~6~7~8]^T$, $\textbf{c}=[9~10~11~12]^T$,
Then $h_1$ is calculated as follows $(a_1+b_1)c_1=(1+5)\times 9$,
$h_2 =(a_2+b_2)c_2=(2+6)\times 10$, $h_3 =(a_3+b_3)c_3=(3+7)\times 11$, $h_4 =(a_4+b_4)c_4=(4+8)\times 12$.

Finally $\textbf{h}=[h_1, h_2, h_3, h_4]^T$.

Question:
1. $\textrm{E}\left[\textbf{h}\textbf{h}^H\right]$=? in terms of $\textrm{E}\left[\textbf{c}\textbf{c}^H\right]$

2.Can I write $\textbf{h}$ in terms of $\textbf{a}$ $\textbf{b}$, $\textbf{c}$?

You could define a diagonal matrix $\mathbf{D}$ with diagonal elements $d_{ii}=a_i+b_i$. Then vector $\mathbf{h}$ can be written as

$$\mathbf{h}=\mathbf{D}\mathbf{c}$$

and the expectation is

$$\mathbb{E}[\mathbf{h}\mathbf{h}^H]=\mathbf{D}\mathbb{E}[\mathbf{c}\mathbf{c}^H]\mathbf{D}^H$$

(I assume that $\mathbf{a}$ and $\mathbf{b}$ are deterministic and $\mathbf{c}$ is random, because you said you wanted the expectation in terms of $\mathbf{c}$.)


What is the relationship between scales and frequency in a Morlet wavelet transform? I'm using PyWavelets, with a complex Morlet wavelet. Its complex wavelet transform function requires scales as one of its parameters, rather than frequencies. However, I don't understand the relationship between the two. I'm generating my scales like this:

base_scale = CENTER_FREQ * quality / base_freq
far_scale = base_scale / 2**num_octaves
scales = numpy.geomspace(base_scale, far_scale, num=num_octaves*voices_per_octave+1, endpoint=True)

That is, if I expect my signal to be at about 20 kHz, I want to set the scales to be such that the frequencies returned are between 10 kHz and 40 kHz. In this case, the num_octaves variable might be 2, base_freq would be 10000, and voices_per_octave might be 50. Then scales would be an array of length 101, ranging logarithmically from the scale corresponding to 10 kHz up to the scale corresponding to 40 kHz.

However, this only works if CENTER_FREQ (the wavelet's central frequency) is a small number. If I increase it beyond 31, the frequencies returned by:

coeffs, freqs = pywt.cwt(data, scales, wavelet, 1/quality)

(where quality is the sampling rate of the data) are dead wrong. Is there any easier way to determine the needed scales that correspond to the frequencies I want?

I found some details, for those who might come after me.

Scale and frequency should be inversely proportional, which is as one might figure. However, there seems to be an issue (possibly a bug) in PyWavelets that makes this not always the case.

There is an internal function in the PyWavelets library called scale2frequency() that takes 2 arguments plus an optional 3rd: wavelet, scale, and then precision. By its name, I assume that it's made to take a scale and return the associated frequency. However it's an involutory function: $f(f(x)) = x$. So it can also be used to calculate scale from frequency. (This might be undocumented behavior, and I have no idea if it'll still work this way in future versions.)

The problem I was running into involved the optional precision argument. It would appear that as the wavelet's central frequency increases, the required precision became very high. I was able to mitigate the problem somewhat by increasing precision beyond the default value of 8. But as it gets higher and higher, the amount of memory and processing power increases drastically, to the point where I burned through all my system's virtual memory long before it was high enough to return valid answers.

My final workaround involved altering my data's sampling rate by a factor of 1000 or so. In essence, I just told it that my data was in the Hz range rather than the kHz range. This made it so that the code wasn't returning any obvious errors. It's not ideal, but it seems to work well enough.


Question about a transient response statement (Oppenheim-Schafer"s discrete-time signal processing book, 3rd ed) This is the statement from the book:
Eq. (129) only ensures that $y_t[n]$  will be bounded  for all values of $n$.
For example, this may include a case where $y_t[n]$ oscillates
and hence never becomes small. Hence, it is not clear how one can conclude that the transient response must become increasingly smaller as $n\to\infty$ when Eq. (129) is satisfied.

It's true that Eq. $(129)$ just presents an upper bound, from which it can't directly be concluded that the transient response $y_t[n]$ decays to zero. However, stability requires that
$$\sum_{k=0}^{\infty}|h[k]|<\infty\tag{1}$$
which can only be true if $\lim_{n\to\infty}h[n]=0$, i.e., if the impulse response decays to zero. If this is the case, then the transient response
$$y_t[n]=-\sum_{k=n+1}^{\infty}h[k]e^{j\omega (n-k)}$$
must also decay to zero for $n\to\infty$.


How can I remove shadows from an image? I have this image

I would like to remove the shadow from the image. I know a lot of different methods like certain morphological operations have been used to remove shadows:

I have created this mask for the same image

Are there some other methods I could try using this mask that I have created?

EDIT:

input image and mask of same size as requested :

EDIT 2: i generated a 1D invariant image however its not perfect

  I = imread('shadow.jpg');
       J = im2double(I);

      R = J(:,:,1);
      G = J(:,:,2);
      B = J(:,:,3);

     [len,wid] = size(R);

     % Generation of 2-D Log Chromaticity Image.
     for i = 1:len
        for j = 1:wid
           if ((R(i,j)*G(i,j)*B(i,j))~= 0)
              c1(i,j) = R(i,j)/((R(i,j)*G(i,j)*B(i,j))^(1/3));
              c2(i,j) = G(i,j)/((R(i,j)*G(i,j)*B(i,j))^(1/3));
              c3(i,j) = B(i,j)/((R(i,j)*G(i,j)*B(i,j))^(1/3));
           else
              c1(i,j) = 1;
              c2(i,j) = 1;
              c3(i,j) = 1;
        end
    end
end

rho1 = mat2gray(log(c1));
rho2 = mat2gray(log(c2));
rho3 = mat2gray(log(c3));

X1 = mat2gray(rho1*1/(sqrt(2)) - rho2*1/(sqrt(2)));                                         %(1/sqrt(2); -1/sqrt(2); 0)
X2 = mat2gray(rho1*1/(sqrt(6)) + rho2*1/(sqrt(6)) - rho3*2/(sqrt(6)));   %(1/sqrt(6); 1/sqrt(6); -2/sqrt(6))

theta = 120;

InvariantImage = cos(theta*pi/180)*X1 + sin(theta*pi/180)*X2;
imagesc(InvariantImage); colormap(gray)

Cannot understand what i am doing wrong out here please help?

There are dozens of publications dealing with shadow-detection, generating shadow-masks, and indeed some that actually remove shadows - such as the ones mentioned in the previous posts. I could add some to the list, if needed. The problem, however, is IMHO far from being solved. For a quick-start, given a shadow-mask, I suggest (and have tried in the past) the two following approaches. They definitely reduce shadows - just not always seamless, and I'm sure there are publications (not by me) dealing with shadow-removal in similar ways.

Gradient Domain Manipulation Techniques as described here (C and Matlab Codes provided):
http://www.umiacs.umd.edu/~aagrawal/ICCV2007Course/index.html
The gradient-integration approach can be used for a number of image-processing problems, see the slides/presentation for further examples.

The general idea:

Compute the spatial derivatives (gradient-images) for all color-channels.
Use the shadow-boundaries from the shadow-mask to generate a weight-mask that is close to zero on shadow-boundaries and increases to one within a specified neighborhood along the shadow-edge, i.e. orthogonal to a given edge-point.
Multiply the weight-mask from (2.) with all the gradient-images to reduce/damp the shadow-boundaries/edges.
Integrate the gradient-images using the code from the above link.
For RGB-images, from my experience, compute the average of the separate channels of the original images, and scale the integrated images to match these values to prevent "funny" color-artifacts.

Brightness manipulation in the original image-domain.

Use the shadow-mask to generate a weight-mask that is one outside shadow regions, has a smooth transition (upwards) across the shadow-boundary and a larger than one scale-factor inside shadow-regions. As suggested in a previous post, the scale-factor can be estimated from the area immediately surrounding a shadow-area using the average-brightness together with the average-brightness of the shadow-region.
Multiply the original images (per channel) with the weight-mask, subject to clipping.

I have also tried using different color-models, e.g. HSV, that present luminance or brightness directly, which can then be modified independent of color (hue/saturation). This works essentially like the brightness-manipulation, i.e. generate a smooth weight-mask and multiply it with the luminance-channel. Maybe the two approaches, i.e. gradient-integration and brightness-manipulation, can be combined in a clever way, but someone has probably tried that before too.

Hope that helps, kind regards, Derik.


Why does circular symmetry of complex envelopes mean that $\mathbb{E} \Big[\mathbf{x}[i] \mathbf{x}^T[i] \Big] = \mathbf{0}$? I am making my way through this paper here, and I understand everything up to equation 16, and I am stuck here.

Long story short, the vector $\mathbf{x}[i]$ here is an $N$ by 1 vector, consisting of the complex envelopes of signal + noise + interference across $N$ different sensors, taken at some time index $i$. We have a desired signal, an interference signal, and AWGN noise at each sensor. Naturally, the desired and interference signals at each sensor are delayed versions of each other, but the AWGN at each sensor is independent of the AWGN at another sensor.

What I do not understand, is the author's claim, where he states that since the desired and/or interference signals are just delayed copies of each other, "the circular symmetry of the complex envelopes thus implies that:"

$$
\mathbb{E} \Big[\mathbf{x}[i] \mathbf{x}^T[i]  \Big] = \mathbf{0}
$$

Circular symmetry of the complex envelopes? What does that mean, and why do we get the above equation from that?

A complex random variable with dimension being 2 $Z=X+jY$ can be defined as a Random vector with components $X$ and $Y$ which take values along 2 different dimensions. The Probability Density Function of $Z$ is defined using the PDFs of $X$ and $Y$. The Expectation is defined as $\mathrm{E}(Z)=\mathrm{E}(X)+j\mathrm{E}(Y)$ and the variance is defined the usual way as $Z$ being a complex function of $X$ and $Y$ and the Variance is $\mathrm{VAR}(Z)=\mathrm{E}(|Z^2|)-|\mathrm{E}(Z)|^2$. Intuitively speaking, $Z$ is a complex function of two random variables which just assumes value $X+jY, \forall x,\forall y$. The PDF is written as follows: Write the 2-dimensional Random variable as $$\tilde{Z}=\left[\begin{array}{ccc}X \\Y\end{array}\right]$$where $\tilde{Z}$ is a Real random vector. The PDF of $Z$ is the PDF of $\tilde{Z}$.
Now with Complex Random vectors, $\mathbf{Z}=\mathbf{X}+j\mathbf{Y}$ is a vector of Complex Random variables which has a lot of interesting properties in a sense that it has more to offer than a Vector of Real Random Variables. Particularly of importance are the Covariance, Pseudocovariance(this is something new) and the Expectation. In our situation, we have a special Complex Random vector called the Symmetric  Random Vector $\mathbf{Z}$.
Define the $2n$-dimensional Complex Random Vector as $$\tilde{\mathbf{Z}}=\left[\begin{array}{ccc}\mathbf{X} \\\mathbf{Y}\end{array}\right]$$ where $\mathbf{X}$ and $\mathbf{Y}$ are $n$-dimensional real random vectors. The PDF definition follows the same from Complex Random Variable. You can look up what the Pseudocovariance and Variance mean for a Complex Random Vector. I'm referring to Digital Communications by Proakis here in my definitions and notations.

Now as in the comments, if $\mathbf{Z}$ is a Symmetric i.e the PDF doesn't change on rotation, we can say that $\mathrm{E}(\mathbf{ZZ^\mathrm{T}})=0$. This is a consequence of the following. Since $\mathbf{Z}$ and $\mathbf{Z}e^{j\theta}$ have the same PDF, $$\mathrm{E}(\mathbf{Z})=\mathrm{E}(e^{j\theta}\mathbf{Z})=e^{j\theta}\mathrm{E}(\mathbf{Z}), \forall \theta$$ If we had $\theta=\pi$ (the condition has to be satisfied for all $\pi$), we can see that $$\mathrm{E}(\mathbf{Z})=\mathrm{E}(-\mathbf{Z})=-\mathrm{E}(\mathbf{Z})$$ which clearly says that $\mathrm{E}(\mathbf{Z})=0$. Going further, we can prove that $\mathrm{E}(\mathbf{Z}\mathbf{Z}^\mathrm{T})=0$ the same way assuming a special case of $\theta=\frac{\pi}{2}$. And hence Circularly Symmetric Complex Random Vectors are Zero-mean and Proper. The last equation left unproven here is the definition of being proper.


Delay Search of two signals I need to find the timedelay between two signals.  I guess MSE is the best approach. How to find that MSE?  I am writing the code and I have access to the autocorrelation matrix of one signal and the crosscorrelation of both signals plus the power of the signal.  Is this enough information to calculate the MSE?  Any tips on how that would be done?

Thanks.

As long as the two signals are simply delayed versions of each  other the delay is simply given by the point in time where maximum of the cross correlation occurs.

If the signals are uncorrelated, then the question is meaningless.

The tricky part is when the two signals are partially correlated or are filtered versions of the same original signal. In this case you will first need to define what exactly you mean by "delay". One possible way to approach this would be to time shift the signals against each other, subtract them and sum the squares. The shift that results in the minimum energy in the difference would be considered the delay. Something like $d = min[sum((x(n)-y(n-d))^{2})]$


System Properties; Linear, Causal, Time-Invariant, Stable? I know the answers of the below question but I dont know why, except linear, could you explain the rest? Small tick: correct, small tick with dash incorrect answer.

Such questions have been asked here before, and the concepts are simple and important. Dilip Sarwate is right to ask you to add details as to where you get stuck, so we can give better answers (instead of just solving the problem for you).

I'll try to explain the basics in a clear way so that this answer may also be useful for other users, and to help you understand and solve the problem yourself. Assuming you understand linearity, let's move on to the other properties:

time-invariance: if $y[n]$ is the response to input $x[n]$, $y[n-m]$ must be the response to $x[n-m]$ for any (integer) $m$. By inspection, this is obviously the case for systems A, B, and D. System C can get slightly confusing for beginners, but just use the criterion to see that for $x[n-m]$ the output is $x[-n-m]$ which does not equal $y[n-m]=x[-(n-m)]=x[-n+m]$. So system C is NOT time-invariant. For system E a simple substitution of the summation index shows you that the system is indeed time-invariant.
causality: this is actually very simple. Just answer the question "does the output signal at any time depend on future values of the input signal?" If the answer is no, then the system is causal, otherwise it isn't. Systems A, B, and D only use the current input value to compute the output, so consequently they are all causal. System E uses the current and past input values, so it's also causal. Only system C is not causal because for negative values of $n$ the output depends on future values of the input.
stability: usually the concept of bounded-input-bounded-output (BIBO) stability is used, which means that if the input signal is bounded, i.e. $|x[n]|\le M$, then also the output signal is bounded, i.e. $|y[n]|\le N$, for some arbitrary positive constants $M$ and $N$. For systems A, B, C, and D this is obviously the case, because the output can never be unbounded for bounded $x[n]$ (assuming that the constants $A$ and $B$ are finite, of course). For system E that's different, because you have an infinite sum which needs to converge for the output to be finite. But for a bounded input signal you can always construct a case where the output is unbounded. Just choose a constant input signal with any finite value, and the sum will grow without bounds. Hence, system E is not BIBO stable.


Relation between the DTFT and the spectrum of a sampled signal In the $\rm DTFT$ (Discrete Time Fourier Transform) the spectrum is periodic with period of $2\pi$ . A continuous signal when sampled has a spectrum which is a repeated version of its original spectrum before sampling with a period of sampling frequency. Both have periodic spectrum and both are discrete in nature.When we sample continuous signal, its values at 'nTs' is multiplied by impulse fn.But in discrete aperiodic signal it is of finite height( where impulse has infinite height)

Question is:

Can both signals be considered similar because their spectra are having a kind of similarity?
Is there any application for the above observation?

If you have a continuous-time signal $x(t)$, then the two signals you're talking about are

$$\begin{align}
x_c(t) &=x(t)\cdot\sum_{n=-\infty}^{\infty}\delta(t-nT) \\
 &=\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT) \\
 &=\sum_{n=-\infty}^{\infty}x(nT)\delta(t-nT) \\ \tag{1}
\end{align}$$

and you define

$$x_d[n] \triangleq x(nT)\tag{2}$$

The first signal given by $(1)$ is technically a continuous-time signal, even though it is only non-zero at discrete times $t=nT$. The reason why it is considered a continuous-time signal is because it can and must be transformed using the continuous-time Fourier transform (CTFT). So $(1)$ is the continuous-time representation of a sampled signal. Eq. $(2)$ is the discrete-time representation of the same signal. Here the sampled signal is represented as a sequence of numbers. You can't apply the CTFT to $(2)$, but you must use the discrete-time Fourier transform (DTFT).

The nice thing is now that the CTFT of $x_c(t)$ given by $(1)$ and the DTFT of $x_d[n]$ given by $(2)$ are identical. So if the CTFT

$$\begin{align}
X_c(j\Omega) &= \int_{-\infty}^{\infty}x_c(t)e^{-j\Omega t}dt \\
   &= \int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}x(nT)\delta(t-nT) e^{-j\Omega t}dt \\
   &= \sum_{n=-\infty}^{\infty}x(nT) \int_{-\infty}^{\infty}\delta(t-nT) e^{-j\Omega t}dt \\
   &= \sum_{n=-\infty}^{\infty}x(nT)  e^{-j\Omega nT} \\
   &= \sum_{n=-\infty}^{\infty}x_d[n]  e^{-j n(\Omega T)} \\
\end{align}$$

and the DTFT:

$$X_d(e^{j\omega})=\sum_{n=-\infty}^{\infty}x_d[n]e^{-jn\omega}$$

we have:

$$X_d(e^{j\omega})=X_c\left(\tfrac{j\omega}{T}\right)\tag{3}$$

In sum, the signals $(1)$ and $(2)$ are just two different representations of the same signal, and their spectra (one defined by the CTFT, the other defined by the DTFT) are identical.


Sum of Squared Intensity Difference for block matching with ORB features, c++ I am trying to implement feature matching algorithm between image reference and image current.

The pipeline is as follows:

Gray scale Images (intensities = [0-255]) are imported.
Extract ORB keypoints in both reference and current images:

Ptr<FeatureDetector> detector = ORB::create();

Obtain blocks surrounding all keyPoints within specified window_size in both reference and current images. NOTE: If the blocks exceed the image upper or lower bounds, the block intensities would obviously be zero (Example: KP_A in image reference).

4.Instead of using existing feature matching algorithms in opencv, I am trying to utilize sum of squared of intensity differences (SSD) in the blocks acquired across ORB keypoints in reference and current images.

Extracting correct features demands implementing crossCheckedMatching() to ensure features are chosen correctly.
Visualize matched features.

Complete code: here.

Problem 1: My final result looks something like this:

Whereas, I have read in few lecture notes and articles, SSD and block matching incorporated with cross checked matching method should provide decent matches. Does these matched features look any good in someone's eyes or awful after cross check matching procedure?

Problem 2: I do not know how to find the correct threshold for evaluating SSD values between two images (Empirical tuning vs. Existing logical approaches?).

Final goal, adapted from the textbook:

use these matched features in RANSAC to filter outliers.
Find Homography (Theory says, 4 pairs of matched points needed!)
Camera Pose Estimation
3D reconstruction of environment

Apologies for such a long question, but I tried to squeeze it as much as I could.

Your results look reasonable. However, you're likely to get a somewhat unreliable result: the good matches are located in a small portion of the image, while a large spread would help in getting overall more accurate estimates (but that's probably due to your images).
Experimenting will work here! Otherwise, you can compute the expectation of the SSD between 2 blocks that are a rightful match under the noise level of your images.

Eventually, bad matches should be eliminated by the RANSAC part, but you may have to do a trade-off between getting many but sometimes wrong matches (loose SSD threshold) or stricter ones (at the risk of not having enough matches o estimate the pose).


Understading sampling theorem and aliasing A real-valued analog signal with a flat spectrum between f = 0 Hz, and fmax is sampled at a sampling frequency fs = 24 kHz. The sampled signal is then processed with an ideal lowpass filter with cutoff frequency = 0.25pi.

How can we specify the range for omega in which aliasing can be torelated and also determine the maximum signal frequency fmax of the analog signal such that no aliasing components are present in the filter output?

An ideal low pass filter has perfect attenuation above the cut off frequency.  So any aliasing distortion in that band will be filtered out.  In your example, this means any content above 3kHz.

For this application fmax can be made to extend up to the sample rate less the cutoff frequency of the low pass, or 21kHz.


I'm confused about the notion of frequency when it comes to arbitrary data In signal processing class we keep learning about how it's so useful to bring things into the frequency domain instead of the time or space domain but I'm very confused about what the motivation is or what the frequency actually corresponds to when it comes to arbitrary data.

Let's say I have some plot of sin(1:100) in MATLAB

and then plot(1:100, fft(sin(1:100)))

What useful information do I get out of the FFT of this function or in general, any function? And for some arbitrary data set like this (not corresponding to sound, pixels, etc.) how should we interpret the frequency?

Eventhough I cannot really get the exact point of your confusion about the concept of frequency and its use in engineering, let me describe it like this.

Space, time and mass are fundamental constructs that we use to measure and describe physical events and objects using their length, volume, speed, duration and position etc...

When you want to describe a rotating wheel, then time duration of a single rotation is termed as a period; $T_0$ in seconds, as the event is periodic in nature, i.e., it repeats in time, ideally forever if not disturbed.

While the period of the rotating event is a sufficient descriptor of it, another very useful and equivalent descriptor is given by the number of rotations per unit time (such as a second), which is then called the frequency, $f_0$ Hertz (Hz.), of the event, which indicates and quantifies its repetition rate.

The usual relation between the period in time and frequency in Hz. is $$ f_0 = \frac{1}{T_0} $$

Period and frequency of spatial events are also described similarly but using meters for the unit of period and cycles per meter (cpm) for the frequency instead. Sub-units such as cycles per centimeter, per millimeter etc. can also be used. Period of a spatial periodicity can also denoted as its wavelength eventhough it's more suitable for traveling or standing waves at a velocity v.

Note, fundamentaly, that eventhough time and space are physical objects (constructs), the concept of frequency is a purely mathematical one.  The importance of the frequency becomes more clear as modern mathematical inventions such as Fourier's theorem and later linear system thoery become ubiquitous devices of engineering and science.

First came Fourier's theorem (Fourier series) which stated that any periodic continous function (signal) can be expressed as a weighted sum of periodic and continuous sinusoids. The theorem extended to include those non-continuous periodic signals (such as periodic pulse trains) as well. Fourier's theorem is later mathematically modified by Dirichlet conditions to be consistent with the rest of the body. Then It's further generalized into continuous and discrete Fourier transforms so that it coulds be used to analyse even more practical signals which are not even periodic.

The main concept is that a periodic signal of fundamental frequency $f_0$ can also be expressed as a superposition sum of harmonic family of sine waves. So this provides the tool necessary to analyse signals by frequencies. What comes next is to analyse systems by frequencies?

And that's achieved by linear system theory, which states that linear (and time invariant) systems posses eigen-functions of periodic sinusoidal waves (complex exponentials more specifically) of arbitrary frequencies. Eigenfunction property of linear (time-invariant) systems make it very convenient to mathematically and practically analyse their behaviour by using periodic input excitations.

So, together with the Fourier's periodic decomposition and linear system theory based eigen-function properties of systems, engineers have a very fundamental method of analysis and design of practical systems.

That explains why concept of freuqency and its methematical extentions to CFS, CTFT, DTFT, DFS, DFT (FFT) are so important, usefull and versatile for modern engineering work.

That being said, nonlinear systems and transient system analysis is also as important (but more difficult) and sometimes even more important than the steady state, frequency based analysis of linear systems.


Determining Fourier Series coefficient for Discrete time I am trying to solve the proof for determining the Fourier series representation of a periodic signal.
I understand fourier series equation for Discrete time which is $$x[n] = \sum^{}_{k=<N>} a_k  e^{jk\omega_0n}$$
In order to find the Fourier series coefficient we multiply both side by $e^{-jr(2\pi/N)n}$ and summation over N terms, That is:
$$\sum^{}_{n=<N>} x[n] e^{-jr(2\pi/N)n} = \sum^{}_{n=<N>} \sum^{}_{k=<N>} a_k  e^{j(k-r)(2\pi/N)n}$$

After that i couldn't understand how can we get the final term:
$$ a_r = \frac{1}{N} \sum^{}_{n=<N>} x[n]e^{jr(2\pi/N)n}$$

Please, anyone, help me to understand this

In the formula displayed for $x[n]$, the exponentials $e^{jk\omega_0 n}$ should be $e^{-jk\omega_0 n}$. I'm assuming that the sum $\sum\limits_{n = \langle N\rangle}$ is extended over all integers $n$ from $0$ to $N-1$. Fix $r\in \{0,\ldots, N-1\}$. Then

$$\sum_{n = \langle N\rangle}x[n]e^{jr\omega_0 n} = \sum_{n = \langle N\rangle}\sum_{k = \langle N\rangle}a_ke^{-jk\omega_0 n}e^{jr\omega_0 n} = \sum_{k = \langle N\rangle} a_k \sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n}.\tag{*}$$

If $k = r$, then each of the terms $e^{j(r-k)\omega_0 n}$ equals $1$, in which case

$$\sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n} = N.$$

If $k \neq r$, let $z = e^{j(r-k)\omega_0}$. Then $z\neq 1$ and $z^N = e^{2\pi j(r-k)n} = 1$. Thus

$$\sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n} = \sum_{n = \langle N\rangle} z^n = \frac{1 - z^N}{1 - z} = 0.$$

In summary,

$$\sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n} = N\delta_{k,r}.$$

Therefore,

$$\sum_{k = \langle N\rangle} a_k \sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n} = \sum_{k = \langle N \rangle} a_k (N\delta_{k,r}) = N\sum_{k = \langle N\rangle} a_k \delta_{k,r} = Na_r.\tag{**}$$

Combining (*) and (**),

$$Na_r = \sum_{n = \langle N\rangle} x[n]e^{jr\omega_0 n} = \sum_{n = \langle N\rangle}x[n]e^{jr(2\pi/N)n}.$$

Hence

$$a_r = \frac{1}{N}\sum_{n = \langle N\rangle} x[n]e^{jr(2\pi/N)n}.$$


Is zero-padding necessary for speech processing to satisfy linear convolution property? In the FFT domain, for example, if I multiply a noise reduction gain to the noisy signal, $ Y(w)H(w) $, this translates to a circular convolution in the time domain of the same length $ y(t) \circledast h(t) $. If I want to convert this to a linear convolution, I should do zero-padding such that the noisy signal y(t) is now twice the original length minus 1 (assuming the noisy signal $ y(t) $ and filtering signal $ h(t) $ are of the same length).

My question is, is obtaining a linear convolution necessary for speech processing? What happens if I just stay with the circular convolution (i.e. NFFT length = signal length)? What if I just want to maintain the same signal length before and after processing? I am quite confused about the importance of satisfying linear convolution property in the context of STFT.

My answer is not from speech processing angle but from a generic signal processing involving Linear and Circular Convolution.

  Is zero-padding necessary for speech processing to satisfy linear
  convolution property?

If the FFT size is $N$, and the length of result of linear convolution size is $\gt N$, then you need to do zero-padding. Otherwise it will result in time-aliasing in time domain. The effect of this is, at the output of FFT you will get an under-sampled representation of the theoretical continuous time DFT signal.
$$
(Y.*H)(e^{j2\pi k/n}) = (Y.*H)(e^{j\omega})|_{\omega = 2\pi k/N}
$$
but
$$
(Y.*H)(e^{j\omega})|_{\omega = 2\pi k/N} \ne (\sum_0^{N-1}p[n]e^{-j\omega n})|_{\omega = 2\pi k/N}
$$
because of time aliasing. Here $p$ is the sequence obtained by linear convolution.

  What happens if I just stay with the circular convolution (i.e. NFFT
  length = signal length)?

Since your FFT length $N$ is equal to signal length $N$, the above time aliasing will happen.

  What if I just want to maintain the same signal length before and
  after processing?

Your FFT size has to be atleast equal to the resultant linear convolved signal length ($2N-1)$. Because that will be the length of signal you will be working with after the point-wise multiplication. So your FFT size has to be at least $2N-1$ before hand. See this question and its answers too Linear and Circular Convolution in Fourier Domain (DFT)


Is $y[n]=x[n] * x[n^2]$ invertible? Is the following system invertible or not?

$$y[n]=x[n] * x[n^2]$$ where $*$ stands for the aperiodic convolution operator. I have not been able to find a mathematically sufficient argument for it...

The system is nonlinear (bilinear in $x$), with a nonlinear law (square) on indices. Odds are the system is not invertible. One can try to prove it in its full generally, or try to find a counterexample. Let us try the lazy way, using properties at hand.

The system is bilinear. Hence, the output for $k x[n]$, $k\in \mathbb{R} $, will be  $k^2 y[n]$. Resultingly, since $(-k)^2 = k^2$,  $k x[n]$ and  $-k x[n]$ have the same output. If they are nonzero, they differ and yield the same result. Verdict: non-invertible.

This is a first way. Another is possible, along the same line, using the indices $n$ and $n^2$. Fat32 just did it, using prime delays. With more generality: take some signal $x[n]$, assuming that:

there exists  at least one non-square index $n_+$, for which $x[n_+] \neq 0$: for instance $x[-2]  =1$ (many negative numbers are not perfect squares),
$x$ is zero at all square indices:  $x[0] = x[1] = x[4] = x[9]=\cdots = 0$.

Then $x[n^2]$ is identically zero, while  $x[n]$ is not and the convolution is zero.


What is difference between outputs of Fourier transform and Fourier series of a periodic square waveform? We can use Fourier transform of an aperiodic signal and Fourier series of periodic signal. But we can use Fourier transform formula for periodic function also.

Now, let us consider a periodic square wave with fundamental period $T$ . Then I want to ask is

  What is difference between the outputs of Fourier transform and Fourier series of a periodic square waveform?

The Fourier transform $X(\omega)$ of a $T$-periodic function $x(t)$

$$ x(t+T) = x(t) \quad \quad \forall t $$

having complex Fourier coefficients $c_n$

$$ x(t) = \sum_{n=-\infty}^{\infty} c_n e^{j 2 \pi n t/T} $$

$$c_n = \frac{1}{T} \int_{0}^T x(t) e^{-j 2 \pi n t/T} \ dt \tag{1}$$

can be expressed as a weighted sum of Dirac impulses, where the weights are given by the complex Fourier coefficients:

$$X(\omega) \triangleq \mathcal{F}\left\{ x(t) \right\} = 2 \pi \sum_{n=-\infty}^{\infty} c_n \delta\left(\omega - \frac{2\pi n}{T}\right) \tag{2}$$


interview questions for DSP engineer I am preparing an interview which may require some knowledges on DSP. Is there any questions, fundamental concepts, and theories that I would not want miss? Thanks a lot for the suggestions.

When I interview people, my questions depend a lot on the CV they've submitted.

If they're fresh from university, I am much more interested in whether they have ever heard of some of the basic image processing concepts: What is their definition of noise and background? What does a Gaussian filter do? And of course, I want to know if they know how to program as opposed to hacking together a last-minute solution for homework.

If they claim to be more experienced, I want to know how well they know what they're talking about. What's their opinion on supervised vs. unsupervised segmentation? What are good tracking algorithms, and why? What informs their choice of a filtering algorithm? How do they deal with stochastic data? And of course, do they know how to program?

Finally, I'm always happy to find out that they did some research into the topic ahead of time, even if they just read my papers and totally misunderstood them. I can always teach them knowledge, but it is hard to teach effort.


Classification from a feature vector I'm quite new to this; I'm try to classify textures as defective or non-defective. I've used a Gabor filter bank with Matlab which outputs a column vector of the Gabor features of an image. I have a data set of non-defective images and defective images.

My question is, what can I now do with this (or these) feature vectors to classify the texture? I've read about many types of classification (SVMs, Neural networts, nearest neighbour etc.), but couldn't find any similar types of implementation to help me get an idea of what I'm doing - perhaps someone could point me towards a practical example to help my understanding. Many thanks.

You have pairs of features-target's which is the first step.

I'm no expert in texture classification, there are probably specialized and/or well-known good methods.

To make it easy for a beginner, I would take a look at the examples for the Python library scikit-learn:
http://scikit-learn.org/stable/auto_examples/index.html#general-examples
and see what looks interesting, for example a few different random forest examples http://scikit-learn.org/stable/auto_examples/ensemble/plot_forest_iris.html

What you need to run these is a matrix with the filter responses features (the same shape as iris.data in the example (number of examples on the first axis and number of features on the second)) and a binary vector of targets (the same shape as iris.target (number of examples on the first axis))
basically you are replacing the "iris" dataset with your own dataset.

to load text data

import numpy as np
features = np.loadtxt('data.csv')

(as a side note: numpy works much like the builtin matrix stuff in matlab)

Try out a few different classifiers with a few different parameters and see what happens with the results http://scikit-learn.org/stable/modules/generated/sklearn.metrics.classification_report.html.

Make sure you split the data into train and test data http://scikit-learn.org/stable/modules/generated/sklearn.cross_validation.train_test_split.html or use cross-validation http://scikit-learn.org/stable/modules/cross_validation.html to avoid overfitting, again there is examples of use under the scikit-learn homepage.