Sampling Rate Effect on Quantizer SNR I am currently reading the 'Coding' chapter on Rabiner & Schafer's Book on speech processing.

In one of its exercises, the reader is given a simple A/D converter using 16-bit uniform quantization. The converter uses a sampling rate of $Fs = 8000 Hz$. Then it asks what is the effect on $SNR_Q$ (quantization signal-to-noise ration) if we double the sampling rate.

From what I understand, doubling the sampling rate has obviously no effect on the quantizer's range $X_{max}$ as well as the quantization step $\Delta$. Furthermore, the theoretical model for $SNR_Q$ on uniform quantizers is built on top of the hypothesis that the quantization error $e(n)$ follows a uniform probability density, so:

$$
SNR_Q = \frac{\mathcal{E}(\sum_n x^2(n))}{\mathcal{E}(\sum_n e^2(n))}
$$

should remain constant, since both its numerator and denominator are calculated off expected values.

Is my thinking correct or does the change of the sampling rate break any of my assumptions?

You are correct. When using uniform sampling, the sampling rate is unrelated to the quantization noise. Only the quantization step and the signal's dynamic range are relevant for calculating the quantization SNR.


When is phase needed in DFT? I observed than people often compute the DFT of a signal using FFT and just observe the peaks in the "modulus VS frequency" diagram (kind of a power spectrum density). Obviously only half of the information appears in this diagram and it is not possible to retrieve the original discrete signal from this plot only.

Question: In practice, in which cases is the "phase VS frequency" plot needed or not?

From the perspective of filtering, a common LTI system will weight each components in the input signal: (1) weighting means changing the magnitude values (it is apparent in multiplication of two complex numbers), thus it is often to observe the magnitudes (this does not mean the phase is not important; (2) of course, phase information is also part of the whole signal, and although the phase information may be also changed, but the phase is not apparent; (3) we can consider a FIR filter, and compare the changes of both magnitude and phase.


Understanding the relationship between the convolutional feature map and the region proposals in Fast R-CNN In the paper Fast R-CNN available here, I am trying to understand the relationship between the region proposals and the convolutional feature map.

So from what I understand, Fast R-CNN takes in an image and a set of region proposals for object detection, with the goal of taking each proposal and giving it a classification from a set of classes, and a bounding box (or none in the case of background regions). More processing can be done after with the region proposals to merge them etc. but in general this is how it works.

What I don't understand is how the region proposals relate to the convolutional feature map. In section 2.2 the authors state that they test their architectures by pre-training networks with some image net architectures, which we will modify by replacing the max pooling layer by a RoI pooling layer, followed by some fully connected sibling layers for regression and classification. My question is what is the relationship between a region proposal, say a region $(x,y, h ,w)$ given by the top left corner and a height and width, and the convolutional feature map, which may have spatial dimensions (and definitely more channels) that differ from the original image, (see the image below)? It would make sense certainly if the feature map had the same spatial dimensions as the original image, but otherwise I'm not sure.

Notation

Lets say you have an input image I with size WxHx3
ROI proposal x0,y0,w0,h0.
When you forward the image through the feature extractor you get a feature map F with size WfxHfXC size.

The relation is a spatial relation
The ROI doesn't affect the number of output channels from the feature map.
The relation is a spatial relation which maps an input ROI to an equivalent spatial patch on the feature maps which has the same number of channels has the feature map. So the ROI in the output feature map will be of size UxVxC the same number of channels the feature map has.
Tracking an input pixel throgh a CNN
Now lets track a the pixel x0,y0 between successive layers of a CNN.
For simplicity i will consider only the following layers types (The same has in VGG 16 which is the feature extractor used in the article, for other layers the same logic can be applied with small variations):

Convolution layer with "SAME" padding, kernel = 3, stride = 1 and denoted by Conv
Pooling layer with some stride = 2 and denoted by Pool

Now what happens when we insert an image to Convolution layer? In other words to which pixel an input pixel x0,y0 is mapped to in the output of the Convolution?
The output size of a convolutional layer is calculated using the following formula: Output size calculation so in our case we will get that the output size is equal to the input size and that the pixel in the convolution response for pixel x0,y0 in the input is mapped to pixel x0,y0 in the output of the convolution.
For a polling layer with stride 2 we will get that if the input size is WxH the output size will be W/2xH/2 so pixel x0,y0 in the input is mapped to floor(x0/2), floor(y0/2)
Summery of tracking (Given our simplified case)

Convolution x0,y0 in the input will be mapped to x0,y0 in the output
Pooling x0,y0 will be mapped to x0/2,y0/2

Converting x0,y0 in image to x0^,y0^ in the output feature map
So if we have a CNN with Conv->Conv->Pool->Conv->Conv->Pool ....
The pixel x0,y0 in the input will be mapped to the pixel floor(x0/2^(num of pool layers)),y0/2^(num of pool layers)
An ROI mapping
Lets take the case in which the feature extractor is VGG16, which the network used in the article.
All the Convolution layers give output size is the same has the input size and the polling size is half the input size.
The input size is 224X224 and the feature map is 7X7 -> the pixel x0,y0 in the input image is mapped to x0/32,y0/32 in the feature map.
The roi x0,y0,w,h is mapped to x0/32,y0/32,h/32,w/32


Implementing an analog approximation of an FIR filter given its impulse response Problem:

I have the impulse response of a matched filter(therefore its phase and magnitude response. See figure below) of a filter, and  I need to implement its response using only off-the-shelf components i.e low-pass filter, delay lines, multiplexers, etc.. I don't require extreme precision but I'd like to approximate its response as exact as possible.

Possible Solutions:

Cascade low-pass and high pass filters to obtain the overall shape of the magnitude response. Use delay lines to implement two taps to give the magnitude response its "sinusoidal" nature as shown by the above figure.
Use a multiplexers/splitters with delays, phase shifters, and attenuators with values according to the filter's frequency and phase response.

I am extremely new to signal processing so I don't know much of what is out there. So please, if anyone could guide me or give me an insight into how to achieve my goal, I'd really appreciate it.

FIR filters typically have a sinc^x type impulse response spectrum or something that looks like it, in some but not all cases. The sinc^x shape is mainly seen when the impulse response h(n), if plotted on cartesian coordinates, look like an even function time shifted where the coefficients are placed symmetrical about the origin. This is a simplistic answer based on a particular type of FIR filter such as a Type I or II. You have not given enough information to give a more specific answer.

So your analog filter would have to look like a sinc^x or some type of comb filter with lobes falling off faster than a typical comb filter. But again this is a simplified answer.

In short it is theoretically impossible to do what you want exactly, or even approximately, but I can't tell you that you cannot somehow come up with a crude approximation, say using delay lines


Fourier transform of certain noisy function So, I have a noisy signal in time domain, $f(t) = t \eta(t)$ where $\eta(t)$ is white noise with variance $\sigma$ and mean zero, and that it has the property $\langle \eta(t)\eta(t') \rangle = \sigma^2 \delta(t-t')$ ($\langle \ldots \rangle$ denotes ensemble averaging). What I want to do is that I want to calculate the square magnitude of the spectrum in frequency domain upon ensemble averaging, $\langle |z(\omega)|^2 \rangle$. To do this, I follow Dave's answer in https://physics.stackexchange.com/questions/53739/magnitude-of-the-fourier-transform-of-white-noise. So, it will be
$$
\langle |z(\omega)|^2 \rangle = \iint e^{i\omega t} e^{-i\omega t'}t't \langle \eta(t)\eta(t') \rangle dt dt'
$$
Using the property stated above, $\langle \eta(t)\eta(t') \rangle = \sigma^2 \delta(t-t')$, I can obtain
$$
\langle |z(\omega)|^2 \rangle = \sigma^2 \int t^2 dt
$$.
Thus, the ensemble averaged square magnitude of my function is independent of $\omega$.

Now I perform the same calculation using MATLAB and I obtained the following pictures.

It seems like when I take ensemble average of the lower curve, it will be dependent on $\omega$, which contradicts my calculation. Can someone please point out my mistake(s)?

clear all;
clc;

n = 2048;
ft = rand(1,n)*(2+2)-2 + 0;
t = linspace(-10,10,n);
Dt = 1/(t(3)-t(1));

fw = ifft(fft(t.*ft));
w = (-length(t)/2:length(t)/2-1)*2*pi*Dt/length(t);
subplot(211)
plot(t,t.*ft)
xlabel('t')
ylabel('f(t)')

subplot(212)
plot(w,abs(fw).^2)
xlabel('\omega')
ylabel('|z(\omega)|^2')

That's my code. By the way I think you are using normal distribution for the noise, I used white noise type.

After thinking more about this, I am pretty sure that what you really want to do, and what your MATLAB probably in fact is doing, is study the a function like

$\displaystyle f(t)=\int^tdt^\prime (t-t^\prime)\eta(t^\prime)$

with $\eta(t)$ having the correlation property you've already mentioned.  This will give an FT that is definitely frequency dependent and, I think, looking like your MATLAB results.


Multiband processing in one pass I want to split signal in multiple bands and process each band. Number of bands is not yet determined, but it would be more then 6 and less then 11 (depending on the collected results)

As demonstrative example, let's consider common linear PCM sound data file input.wav and I want to calculate RMS per each band. Using a tool such as SoX I can extract only one band at a time (1000-2000 Hz in this example):

sox input.wav output.wav gain -n sinc 1000-2000

and pipe output to processor and extract wanted feature.

What I would appreciate is solution to process all previously determined bands in one pass.

It doesn't need to be a ready made program like SoX, but preferably Matlab/Octave script or Python/SciPy/NumPy and similar approach.

The most direct way to make these calculations is to design a set of filters to your specifications, then read the data out of the WAV files and run it through the filters. In Python/SciPy, you can use Python's wave module to read and write WAV files, SciPy's firwin to design a windowed sinc filter like SoX's sinc effect, and SciPy's lfilter to apply each filter. You will have to design each filter separately and call lfilter for each filter.

A simple way to approximate the energy in several bands is to first calculate the Discrete Fourier Transform (DFT) of consecutive frames, and sum up the energy in the bins corresponding to your bands of interest. The DFT approximates a signal with a sum of harmonically related sinusoids. The advantage of this technique is that it's easy to sum up several arbitrarily spaced bands. A big disadvantage is you have to apply a window to control spectral leakage. This makes the energy calculation an approximation. You can use NumPy's rfft function to calculate the Fourier Transforms.

If your bank of filters is related, for instance all the same bandwidth, or centered at harmonically or logarithmically related frequencies, there are some clever techniques that fall under the concept of filterbanks. This is more specialized, and I don't know of much pre-packaged code you could use.


What is the unit of x-axis of a sound signal? I'm confused about the unit of x-axis.

The below signal has a duration of 3 second and sampling frequency 44100. So should i write Time(s) for the unit of x-axis.

However, in the graph, it shows the number of samples i.e 44100 * 3 = 132300.

The below signal has a duration of 3 second and sampling frequency 44100. So should i write Time(s) for the unit of x-axis.

No. In your graph, the X axis seems to be simply the sample number. I mean, it's a bit strange to ask us, considering you generated the graph, but if you benefit from having it spelled out:

In DSP, signals are represented by a sequence of samples. And if you plot these samples, there's no inherent "units" on the X axis. Your signal has been reduced to a sequence of digital numbers, and the fact that there was a sampling rate that implies some period of time between each sample is simply not there in DSP.


Calculating Haar features I am confused on how to solve this problem, does anyone have any ideas?
https://www.merl.com/publications/docs/TR2004-043.pdf
For a Haar filter $H \in \mathbb{R}^{3 \times 5}$ defined as below:
What is the equation to calculate the Haar features using integral image II?
In this problem an input images is $I \in \mathbb{R}^{m \times n}$ and integral image is denoted by $II\in \mathbb{R}^{m \times n}$
$H = \begin{bmatrix}
0 & -1 & -1 & 1 & 0\\
0 & -1 & -1 & 1 &0\\
0 & -1 & -1 & 1&0\\
\end{bmatrix}$

As mentioned the Haar filters use a lookup table with the say $A_{i,j} = \sum_{k=0}^{i} \sum_{l=1}^{j} I_{k,l}$
Now consider a the convolution with a filter $H \in \mathbb{R}^{m \times n}$, such that $H_{i,j}=1$, the convolution will can be calculated as
$$\begin{eqnarray}Y_{i,j} &=& \sum_{k=1}^{m} \sum_{l=1}^{n} H_{k,l} I_{i-k+1, j-k+1} \\
&=&  \sum_{k=1}^{m} \sum_{l=1}^{n} I_{i-k+1, j-k+1} \\
&=&  \sum_{k=i-m+1}^{i} \sum_{l=j-n+1}^{j} I_{k, l} \\
&=&  \sum_{k=i-m+1}^{i} \left( \sum_{l=0}^{j} I_{k, l} -  \sum_{l=0}^{j-n} I_{k, l}\right) \\
&=&  \sum_{k=0}^{i} \left( \sum_{l=0}^{j} I_{k, l} -  \sum_{l=0}^{j-n} I_{k, l}\right) - \sum_{k=0}^{i-m+1} \left( \sum_{l=0}^{j} I_{k, l} -  \sum_{l=0}^{j-n} I_{k, l}\right)\\
&=&  \left(A_{i,j} -  A_{i,j-n}\right) - \left( A_{i-m,j} -  A_{i-m, j-n}\right)\\
\end{eqnarray}
$$
Of course you can generalize this for $H_{i,j} = c$ for constant $c$ (by multiplying the resulting image), and you can use the shifting properties of the convolution if the filter is padded.
Now looking to your $H$, it can be decomposed in two rectangular regions constant regions. It can be decomposed as
$$H_1 + H_2 = \begin{bmatrix} 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 & 0\end{bmatrix} + (-1) \begin{bmatrix} 0 & 1 & 1 & 0 & 0 \\  0 & 1 & 1 & 0 & 0  \\  0 & 1 & 1 & 0 & 0 \end{bmatrix}$$
To convolve with a the non-zero submatrix of $H_1$, a matrix in $\mathbb{R}^{3 \times 1}$ you will use $A_{i,j} + A_{i-3, j-1} - A_{i, j-1} - A_{i-3,j}$
To convolve with the non-zero submatrix of $H_2$, a matrix in $\mathbb{R}^{3 \times 2}$ you will use $A_{i,j} + A_{i-3, j-2} - A_{i, j-2} - A_{i-3,j}$
Of course you can have other conventions for the indices...


Inverse Fourier transform Of a triangular impulse I have to find the expression of this graphic and after find the inverse Fourier transform of it. First of all I found that the expression of the graphic is $$ X(f) = \frac{1}{2} tri (\frac{f+f_0}{B}) - \frac{1}{2} tri(\frac{f-f_0}{B})$$. Now to find inverse Fourier transform , my book give me the advice to multiply numerator and denominator for i. But $$ i=e^{i \frac{\pi}{2}} $$ and $$ -i=e^{-i \frac{\pi}{2}} $$. Now I know that the Fourier transform of a triangular impulse is $$ (sinc(f)^{2}) $$ and that $$ \frac{d}{dt} tri(t) = rect ( t + \frac{1}{2}) - rect ( t - \frac{1}{2}) $$ but I don’t know how to apply correctly integration property of my x(t). Thank you so much

It's nice to see alternative paths, like this one: First of all, properties:

Frequency shifting: $x(t)e^{j2\pi f_0 t} \longleftrightarrow X(f-f_0)$
Duality: if $x(t) \longleftrightarrow X(f)$ then $X(t) \longleftrightarrow x(-f)$

You have correctly expressed the depicted transform: $$X(f) = \frac{1}{2}\mathrm{tri}\Big(\frac{f+f_0}{B}\Big) - \frac{1}{2}\mathrm{tri}\Big(\frac{f-f_0}{B}\Big)$$

As you said, $\mathrm{tri}\Big(\frac{t}{B}\Big) \longleftrightarrow B\mathrm{sinc}^2(fB)$. From the duality property, you have $$B\mathrm{sinc}^2(tB) \longleftrightarrow \mathrm{tri}\Big(\frac{-f}{B}\Big)$$ but since $\mathrm{tri(\cdot )}$ is an even function, you can write $$B\mathrm{sinc}^2(tB) \longleftrightarrow \mathrm{tri}\Big(\frac{f}{B}\Big)$$

Now you have your FT-pair you need. If you apply the frqeuency shifting property on $\mathrm{tri}\Big(\frac{f\pm f_0}{B}\Big)$, you can easily get... what? :-) You can continue from this point.


Which kind of machine learning algorithm is Hidden Markov Model? I have used Hidden Markov Model algorithm for automated speech recognition in a signal processing class. Now going through Machine learning literature i see that  algorithms are classified as "Classification" , "Clustering" or "Regression".
Which bucket does HMM fall into? I did not come across hidden markov models listed in the literature.

I would be tempted to reply "none", or "both classification and clustering".

Why "none"? Because HMMs are not in the same bag as support vector machines or k-means.

Support vector machines or k-means are specifically designed to solve a problem (classification in the first case, clustering in the second), and are indeed just an optimization procedure to maximize an "expected goodness of classification" or "goodness of clustering" criterion. The beauty lies in the choice of the criterion and optimization procedure. HMM are not an algorithm per se. They are a specific kind of probability distribution over sequences of vectors - for which we know good parameter estimation and marginal distribution computation algorithms. But asking whether they are in the "clustering" or "classification" family is as ridiculous as asking whether the Gaussian distribution is supervised or unsupervised learning.

Why "both classification and clustering"? Because of the following: Being probability distributions, HMM can be used for classification in a bayesian framework; and being model with hidden states, some latent clustering of the training data can be recovered from their parameters. More precisely:

HMM can be used for classification. This is a straightforward application of the bayesian classification framework, with the HMM being used as the probabilistic model describing your data. For example, you have a large database of utterances of digits ("one", "two", etc) and want to build a system capable of classifying an unknown utterance. For each class in your training data ("one", "two", you estimate the parameters of a HMM model describing the training sequences in this class - and you end up with 10 models. Then, to perform recognition, you compute the 10 likelihood scores (which indicate how likely the sequence you want to recognize has been generated by the model), and the model with the highest score gives you the digit. In the Rabiner tutorial on HMMs, the training stage is "Problem 3", the classification stage is "Problem 2".

HMM can be used in an unsupervised fashion too, to achieve something akin to clustering. Given a sequence, you can train a $k$-state HMM on it, and at the end of the training process, run the Viterbi algorithm on your sequence to get the most likely state associated with each input vector (or just pull this from the $\gamma$ during the training process). This gives you a clustering of your input sequence into $k$ classes, but unlike what you would have obtained by running your data through k-means, your clustering is homogeneous on the time axis. For example, you can extract the color histograms of each frame of a video sequence, run this process on this sequence, and you'll end up with a break-down of the video into homogeneous temporal segments corresponding to scenes (the unpractical bit is that you have to set up the number of scenes $k$ in advance). This technique is commonly used in automatic, unsupervised, structure analysis of video or music.


What is the right algorithm to detect segmentations of a line chart? To be concrete, given 2D numerical data as is shown as line plots below. There are peaks on a background average movement (with small vibrations). We want to find the values of pairs (x1, x2) if those peaks drops down to average; or (x1) only if the line doesn't back to the average.

There are thousands of such 2D data.

What is the right statistic or machine learning algorithm to find x1 and x2 above without plotting?

Basically, it looks like you need to know some threshold to detect peaks:

detect bottom edge (Xbe)
calculate the moving average (until the last edge)
detect peak itself (current value - average value > threshold)
start of peak found: X1=Xbe
detect peak end (current value - average value < threshold)
detect bottom age (Xbe)
end of peak found: X2=Xbe

Bottom edge detection in the simplest way is when you make moving subtract for two consecutive values and sign changes from negative to positive.


Spectral analysis of positive signals Suppose that I have a sensor that can acquire samples $X[k]$ of the Fourier transform of an unknown signal $Y[t]$. An example is MRI, where the acquired data is in $k-$space. Now suppose that the unknown signal $Y[t]$ is known to be real and non-negative. My question is: is there a principled way to incorporate this knowledge into the spectral analysis algorithm that will estimate $Y[t]$ from $X[k]$, in order to produce an estimate with less bias or variance? I am thinking at non-parametric spectral estimation algorithms. A naive way of course would be to take the real part of $Y[t]$ and clip the negative values, but this does not seem to be optimal. I am looking for some sort of Cadzow's denoising method for spectral data.

To give a complete answer to this question you're going to need to provide more details about the kind of models you're considering in the first place. But yes, in many cases you can augment those models with a priori constraints on $Y[t]$, such as $0 \leq Y[t] \leq 1$.

For example, if the standard model has some sort of least-squares structure, then adding constraints of that type turns the problem into a bound-constrained least squares problem. There are a variety of approaches to solving such problems, and while they are more expensive than standard least squares, they are quite tractable. And it's very likely that such constraints will produce a better reconstruction.

Even without knowing more, though, I will say this: if your modeling approach does not produce real signals naturally, then you are almost certainly using the wrong modeling approach. It concerns me that you are even proposing taking the real part of the output of some other model. You should be searching the space of real signals if you know that to be the underlying structure.


Rate of fall of sidelobe level Anyone can explain me about how measure side lobe roll-off rate of kaiser window from Spectrum in matlab where side lobe roll-off rate is the asymptotic decay rate, in decibels per octave of frequency.

How adjust Frequency axis if length of FFT = "500" by padding "473" zeros  with length of window= "27" to frequency in Hz.

It would be helpful to study the following 2 references:

[1] F. J. Harris, “On the Use of Windows for Harmonic Analysis with the Discrete Fourier
Transform,” Proc. IEEE, vol. 66, no. 1, Jan. 1978, pp. 51-83.  It is also available at:
http://web.mit.edu/xiphmont/Public/windows.pdf

[2] L. R. Rabiner, B. Gold, Theory and Application of Digital Signal Processing, Englewood Cliffs, NJ, Prentice-Hall, 1975. (Section 3.11 on pps. 93-94 refer to the Kaiser window, and Figs. 3.15 (p. 96) and 3.16 (p. 97) show a time and frequency domain 257 point Kaiser window, respectively).

Reference [1] above is the classic windows paper by Harris.  It is a good learning experience to compute the examples of windows from that paper and try to duplicate his graphs.

First order of business – is it a Kaiser window as shown in ref. [2], or a Kaiser-Bessel window, as shown in ref. [1]?  No problem – both authors use Bessel functions when calculating the window, so either term is probably correct, depending on who you talk to.

But wait a minute – the window calculation in [1] (eq. 46a, p. 73) is different than in [2] (eq. 3.54 on p. 94)!  One of them goes from -(N-1)/2 to (N-1)/2 and the other goes from 0 to <= |N/2|.  If the former is to be calculated correctly, N has to be an odd number, and for the latter, N has to be an even number. And yet, for both, adding in the zero point means that BOTH windows will have an odd number of total points!  And as if that weren't bad enough, the equations differ in other ways – [1] uses pi*alpha, while [2] uses beta.  And a denominator in [2] uses an N-1, whereas [1] uses N.  So which one is correct?

The answer lies in: what is the window being used for?  In [2], it is being used for a FIR filter.  As clearly noted in [1, p. 52, two paragraphs beginning with: “Since the DFT essentially considers sequences to be periodic, ...], windows for FIR filters have a different symmetry from those used in a DFT, and the latter have to be shifted by N/2 points.  The author of [1] identifies this as “a major source of window misapplication.”

So you'd better be sure that the window you're calculating is the right one (e.g.: eq. 46a in [1], appropriately right-shifted, and with the proper symmetry).

Another thing you have to be careful about is whether your window point calculations are for an odd or even N.  Whenever you see a -(N-1)/2 to (N-1)/2 range, it's probably odd (try to plug an even N into the indexes, and the end points are fractions).  If you see an N/2 somewhere, it's most likely for an even N.

As I recall, when calculating eq. 46a in [1], I programmed it using -(N-1)/2 to (N-1)/2 and then right shifted the result.  I wasn't quite sure I had the symmetry right, but after experimenting with a few different ways of calculating things, I concluded that I was using enough points such that I got only minor differences in the end results when trying different calculations.

Windows can be confusing at times.  Perhaps the best thing to do is to read [1] again and again and again – there are all kinds of things that you may not understand at first, and it takes time to become comfortable with some difficult concepts.  And it helps to know about papers written by other authors that correct some of the minor errors in [1].

As for zero-padding, if you do so (presumably by appending zeros to your time domain window), it will provide you with an interpolation of your frequency domain points.  I usually do that because it allows me to see the sidelobe structure more clearly.  Padding by 2 or 4 or 8, etc. is typical.

As for plotting – that's another area where you can make a lot of mistakes.  First, compute magnitude squared of your FFT results.  Then, normalize the data such that all your amplitudes are in dB relative to W(0).  Something like 10log10[W(n)/W(0)] (depending on how you're scaling things).  Your highest plot point should be W(0) at zero dB, and all the rest negative.  This is the way windows are almost always depicted in the frequency domain.

As for the frequency axis, well, you can plot versus log(f) – it would make it easier to see the asymptotic decay of the sidelobes.

Of course, you could always cheat (not recommended) by using the approximation for the frequency response shown in eq. 46b on p. 73 in [1].  Or just look at the table for Kaiser-Bessel on p. 55, which shows a -6 dB roll-off rate for the various Kaiser windows.

P.S. - I have no idea who downvoted this question – seems to me to be the type of thing that a lot of practitioners would want to know.


Android: fundamental frequency I want to find the fundamental frequency for human voice in an Android Application. I want to use autocorrelation method. My code is this:

public double calculateFFT(double[] signal)

     {
      final int mNumberOfFFTPoints =1024;

      double[] magnitude = new double[mNumberOfFFTPoints/2];
      DoubleFFT_1D fft = new DoubleFFT_1D(mNumberOfFFTPoints);
      double[] fftData = new double[mNumberOfFFTPoints*2];
      double max_index=-1;
      double max_magnitude=-1;

      final float sampleRate=44100;
      double frequency;

      for (int i=0;i<mNumberOfFFTPoints;i++){

       //fftData[2 * i] = buffer[i+firstSample];
       fftData[2 * i] = signal[i];  //da controllare
       fftData[2 * i + 1] = 0;

       fft.complexForward(fftData);
      }

      for(int i = 0; i < mNumberOfFFTPoints/2; i++){

       magnitude[i]=Math.sqrt((fftData[2*i] * fftData[2*i]) + (fftData[2*i + 1] * fftData[2*i + 1]));

       if (max_magnitude<magnitude[i]){
        max_magnitude=magnitude[i];
        max_index=i;
       }
      }
      return frequency=sampleRate*(double)max_index/(double)mNumberOfFFTPoints;

 }

I used "edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;"
Is my "return"  the value of fundamental frequency or not?

thanks

I changed my code in this mode:

    public double calculateFFT(double[] signal)

    {

        final int mNumberOfFFTPoints =1024;
        double[] magnitude = new double[mNumberOfFFTPoints/2];
        DoubleFFT_1D fft = new DoubleFFT_1D(mNumberOfFFTPoints);
        double[] fftData = new double[mNumberOfFFTPoints*2];

        double max_index=-1;
        double max_magnitude=-1;

        final float sampleRate=44100;
        double frequency;

        for (int i=0;i<mNumberOfFFTPoints;i++){

            fftData[2 * i] = signal[i];
            fftData[2 * i + 1] = 0;

            }

        fft.complexForward(fftData);

        for(int i = 0,j=0; i < mNumberOfFFTPoints/2; i += 2,j++){

            magnitude[j]=Math.sqrt((fftData[2*i] * fftData[2*i]) + (fftData[2*i + 1] * fftData[2*i + 1]));

        }
         fft.complexInverse(fftData, false);

         for(int i=0;i<magnitude.length;i++)

            if (max_magnitude<magnitude[i]){
                max_magnitude=magnitude[i];
                max_index=i;
            }

            return frequency=sampleRate*(double)max_index/(double)mNumberOfFFTPoints;

}

in this mode I calculate autocorrelation thanks inverse transform "fft.complexInverse(fftData, false);"; isn't it?
And also in this case isn't the "return" the fundamental frequency?

Very often not.  Your posted code is for finding an FFT magnitude, not an autocorrelation result, and thus likely to fail at estimating the fundamental pitch of sound/voice/music with typically high overtone or harmonic content.


Measuring similarity between two very similar pictures I know similar questions have already been asked, but as I understand, little differences in my goal make a difference on a decision on which algorithm to use. As I am a total novice in this area, I am asking which algorithm would be best for measuring the similarity of two very similar images - e.g. two photos of faces which were taken from the same angle and are grayscale. Basically, this algorithm would need to quantify the difference or similarity of these two very similar pictures.

What software can I use to do that? Basically, I need something that inputs 2 pictures and then produces a number which tells me how different/similar they are.

EDITED: Picture examples

Hey something as similar as these two pictures  and . And which software can be used for these analyses?

This measure really depends on how your images look like. A very basic thing would be the calculate the point-wise difference between both images and summing them up. The smaller this value is, the more equal they are. Possibly, you can perform some Gaussian Smoothing and contrast equalization before. But, that's really depending on the images.


FFT widow to dsp processor. I am applying 8 samples and triangular window to my fft engine. I am getting the results but not accurate (approximately equal with some tolerance).

I am calculating MATLAB reference as below :

a=[1,2,3,4,5,6,7,8];
b=triang(8);
c=a*b;
d=fft(c)

Is it the correct way to apply window for signal processing unit?

No. You are doing a matrix mutliply and not a vector multiply. Correct version would be

a=[1,2,3,4,5,6,7,8]'; % column vector
b=triang(8); % also column vector
c=a.*b;    % vector multiply
d=fft(c)


audio quality evaluation I am creating a dataset, where in many speakers are recording their speech signals using different microphones and in different environment.  I would like to evaluate the quality of audio and discard the audio of poor quality.  Kindly give leads towards goods measurements (like SNR, PESQ(with reference)), so that a clean dataset can be created.

In specific, I want to check microphone quality from audio and dicard the audio with low quality.

Also please inform about the ideal value of these measurements.

This is the solution I think would work, there might be other more accurate and efficient ways to do this.

You might do this steps for every recording:

Find the silent sections (using the algorithms stated below)
Measure Standard Deviation of those silent portions of your recordings as an indication of noise power.
Select Standard  Deviation of silence sections of one the recordings that you are sure has acceptable quality, as a your threshold value.
Discard any other recording that its estimated noise power is higher that specified threshold (which you estimate in step 3).

For detection of silent portions of the voice, I suggest simply segment your signal into some non-overlapping windows, estimate signal power

$Power=sum(Segment_i^2)/lenght(Segment_i)$

any signal segment which has a signal power below a per-defined threshold, might be considered as silent section.
However this method is not very accurate. If you have enough time, try an algorithm called Voice Activity Detection (VAD), which is more  widely used (the ITU-T G.729 standard uses VAD  to reduce the transmission rate during silence periods of speech). Please, take a look at http://practicalcryptography.com/miscellaneous/machine-learning/voice-activity-detection-vad-tutorial/


Least Squares and Auto and Cross Correlation I am trying to understand why auto and cross correlation helps find the best fit line in least squares. I have an equation as stated here:

$Ax=b$ -- I have not exact solution, so I use the least squares equation.
$$A^TAx=A^Tb$$
I can solve and use the equation but I am trying to figure out why cross and auto correlation helps. Can anyone shed some light on that? For example, what do cross and auto correlation represent?

I am learning about signal processing and linerisation of a signal. $A$ is the input unpredistoted signal and $b$ is feedback distorted signal. $x$ is the coefficients that will be used to weight the Look Up Table. $Ax=B$ is solved in the adaptive algorithm box, in the link below.

http://www.altera.com/end-markets/wireless/advanced-dsp/predistortion/wir-digital-predistortion.html

Cross-correlation and auto-correlation are not directly linked to least squares.

Least squares can create a line of best fit for a set of points.

Cross correlation is a measure of the similarity of one signal to another (as a function of time lag to one signal). The discrete (digital) cross-correlation is an array where each element tells you how much the signals correlate with a certain amount of time lag. This operation is performed naively with a sliding window, or faster using Fourier, or other transforms. xcorr(a,b) = IFFT(FFT(a)* conj(FFT(b)))

Auto correlation is a measure of the similarity of a waveform to itself as a function of time.

There are normalised variants of cross-correlation and auto-correlation. These are created by measuring the standard deviation of each signal and dividing the output by the multiple of both. This enables you to classify the level of correlation much more easily. 0 implies no correlation and 1 implies identical signals.


STFT window width vs. detectable Hz range Assuming 44100 samples per second in a short 4-8 second-long audio .wav file, if I want perform FFT to detect the power spectrum, amplitude, and phase shift from 20 Hz to 10000 Hz (humans lose the ability hear above 10000 Hz with increasing age), I am thinking I will need window lengths of at least 4410 samples?  This is because, within a second, a 20 Hz cycle will occur every 2205 samples (44100/20), but due to Nyquist, I will need twice the 2205 or 4410 samples.

For human hearing, what is the typical window width used so that deconvolution will result in clean audio quality?  For some reason, I think the long wavelengths of 20 Hz will require window widths no shorter than 4410 samples.  Given this, why are their STFT matrices shown in reports which have 1024 bins in the time domain for 4 second .wav files?  (see the image below from the NMF chapter at musicinformationretrieval.com).  The 1024 bins must be based on the jump size -- that site never states that they employ windows widths of 4410 with much shorter jump sizes.

UPDATE

The goal is to replicate the NMF results presented here via OOP with a compilable programming language (.NET) -- not an interpretive language like R, Matlab, Librosa, etc.  Their STFT results were obtained using Librosa by merely invoking the syntax S = librosa.stft(x). It would be futile to chase down Librosa's default settings for FFT bin spacing and window jump length for STFT matrix generation by Librosa because I am not using Librosa.

Regarding the first paragraph of your question: The number of samples used in a discrete Fourier transform (DFT) depends on your desired DFT bin spacing (the frequency difference between adjacent DFT bins) in the frequency domain. If that desired bin spacing is Fo Hz, and the time-domain signal sample rate is Fs samples/second, then the number time samples, N, is:

N = Fs/Fo samples.             (1)

For example, at Fs = 44100, if you select N = 4410 samples then the DFT bin spacing will be:

Fo = Fs/N = 44100/4410 = 10 Hz.        (2)

However, for the traditional radix-2 fast Fourier transform (FFT) N must be an integer power of two. So using the above Eq. (2) when N = 2^12 = 4096 then Fo = 44100/4096 = 10.766 Hz.

Regarding the second paragraph of your question: I don't fully understand that paragraph but if Fs = 44100, someone performing 1024-point FFTs (N = 1024) produces spectra where the FFT bin spacing is:

Fo = Fs/N = 44100/1024 = 43.06 Hz.        (3)

Does that Fo = 43.06 Hz FFT bin spacing seem correct based on other information you see from the NMF chapter?


Butterworth filter's gain formula does not agree with R's `signal` package I'm trying to calculate the Butterworth filter gain. If I use the formula mentioned on Wikipedia:
$$
G^2(\omega) = \frac{G_0^2}
{1+\left(\frac{j\omega}{j\omega_c}\right)^{2n}}
$$
I don't get a matching result from calculating the gain directly from the filter's magnitude using R's signal package.
library(signal)

# Butterworth filter

# Gain formula from wikipedia
Butterworth_gain <- function(freq, cutoff_frequency, n = 1) {
        1/(1+(freq/cutoff_frequency)^(2*n))

}

bf <- butter(n = 2, W = .6, type = "low")
bfr <- freqz(bf)

plot(bfr$f, abs(bfr$h)^2, ty ='l')
lines(bfr$f, Butterworth_gain(bfr$f,  pi*.6, n = 2), col = 'red')

The function butter() computes the coefficients of a discrete-time ("digital") Butterworth filter, whereas the gain formula you used is valid for a continuous-time ("analog") Butterworth filter. According to the R documentation you can use butter() to compute an analog filter using the plane argument.
The gain of a discrete-time Butterworth filter is obtained by the bilinear transform, which substitutes the "analog" frequency variable by $\tan(\omega/2)$, where $\omega$ is the frequency in radians, normalized by the sampling frequency. For a discrete-time unit gain lowpass filter of order $n$, the squared gain is given by
$$G^2(\omega)=\frac{1}{1+\left(\frac{\tan(\omega/2)}{\tan(\omega_c/2)}\right)^{2n}},\qquad |\omega|\le\pi\tag{1}$$


Comparisons of FIR Causal filters of type HRIR in Discrete Time and Frequency Domains I have just implemented a Discrete Time Convolution between HRIR filter at a range of angles. I also implemented an Overlalp-and-Add method and FFT and IFFT to compute the same convolution in Frequency Domain.
The HRIRs are the same measurements which are taken from a public database for both DT and Overlap-And-Add methods.

The input signal is an audio file and the output is supposed to be a 3D effect of producing a sound that moves in the horizontal and vertical planes for defined angle steps for both methods. The horizontal plane movement is a full circle at the level of the listener's ears from left to right and the vertical movement is a half circle from left ear to right ear at the level of the listener.

My questions are:

1) If I wish to compare the filtering done to each HRIR angle above in Discrete Time, is it beneficial to check the shapes, magnitudes of each HRIR angle and compare it against each other angle?

2) With regards to comparing the filtering done in the the Frequency Domain. What parameters should I compare?

3) Does the FFT of the individual HRIR at an angle help to compare the results?

Note that I have only compared the quality of the sound in the output file in both DT and Frequency Domains.

These are the figures of filtering at certain times in Discrete Time and Frequency Domains


Calculate properties of a given Modulated Signal I'm having a hard time trying to understand signal modulation graphics and how can I calculate/find the properties of the signal. The properties I'm talking are:

Max. amplitude
Carrier frequency
the time it takes to code a bit

Can I do it just observing this two graphics? If not, can I do it using matlab?  I'm a little confused since I'm used to see amplitude and frequency graphics but with bars and discrete values, not with continuous functions.

This is one of the examples I was talking and what I believe is a FSK modulation:

Establishing the system's parameters from time/frequency plots with absolute certainty is impossible, but you can make educated guesses.

In your example, the signal seems to be binary FSK; it seems to have continuous phase. Assuming that is correct, then you can try guessing the parameters you're interested in:

Maximum amplitude: from the time plot, you can see that the amplitude is 2. The units are not specified.
Carrier frequency. The carrier could be said to be the frequency at the middle of the two frequency peaks (seems to be around 350 Hz).
Bit interval. It is the time interval for which the frequency stays constant. If we assume that the bits transmitted are 1,0,1,0,1, then it seems like the bit interval is 0.015 seconds or so. Note that the bits transmitted could be 1,1,0,0,1,1,0,0,1,1 with a bit interval of 0.015/2 seconds. This is a case where complete certainty is impossible without more information about the system.


Clarification concerning power spectral density Many books in signal processing,  e.g. Papoulis [1], define power spectral density (PSD) as:

$$S(\omega)=\sum_{k=-\infty}^{\infty}R_{xx}(k)e^{-j\omega k}$$

Which is the fourier transform of the correlation function:
$$R_{xx}(\tau)=E[x(t)x(t-\tau)]$$

However, some authors in the context of time series analysis, for example Jenkins [2] define it as:

$$\Gamma_{xx}(\omega)=\sum_{k=-\infty}^{\infty}\gamma_{xx}(k)e^{-j\omega k}$$

Which is the fourier transform of the covariance function:

$$\gamma_{xx}(\tau)=\hbox{Cov}[x(t)x(t-\tau)]=E[x(t)x(t-\tau)]-E[x(t)]E[x(t-\tau)]$$

As I understand, the Wiener-Khinchin theorem considers $R_{xx}(\tau)$

Can someone clarify why $\Gamma_{xx}(\omega)$ is a valid definition of PSD.

[1] Papoulis, A. (1965). Probability, random variables, and stochastic processes.

[2]  Jenkins,G. Watts,D. (1968), Spectral Analysis and Its Applications

Considering you are using PSD and not BiSpectra implicitly impose you assumed your process is WSS (you only consider the delay between samples and not their absolute time regarding the time origin).

Considering WSS assumption your process' ensemble mean is constant and doesn't change with time (E[x(t)]=Constant). So the only difference between co-variance and auto-correlation is a DC bias which only affect your DC component of your PSD so for zero mean processes(E[x(t)]=0) there won't be any difference between these definitions.

Considering practical estimation of PSD always there is a large DC component in PSD which always leaks to other frequencies and mask them out. So in practice they usually remove the DC component (removes a DC bias from the process then calculate the PSD).


What is the unit of autocorrelation function? In general, for autocorrelation of the deterministic signals,from the formula what is the unit of it.

It depends.  Normalized autocorrelations have no units because the units are divided out as part of the normalization process.

Non-normalized autocorrelations have the original data's units squared.


Convolution of an Image with a Kernel That Is a Product of Two Functions Suppose that $G(i,j)$ is a Gaussian decay function on the distance between points $i$ and $j$ of an image. In addition, $D(i,j)$ is the difference between the VALUES of the image at those points.
Now, at every point $i$ of an image I need to compute the summation of the neighboring $G(i,j) D(i,j) $ , this is, the product of both functions. This is a convolution-like operation: for every point I add the contributions of a function at every neighboring point.
However, I do not know if that can be expressed as a convolution because G and D are functions on different variables: the distance (a raw gaussian filter) and the image intensity difference.
I need help at expressing the summation as a combination of convolution operations in order to easily code it in a program without the need for loops. In addition, it would allow me to possibly do it in the frequency domain.
Any help will be appreciated!

It is equivalent to a bilateral filter. You can move the kernel to each local part of image to do the dot product, then sum them up to get the result in a loop. There is also a fast implementation with O(1) in matlab programming for your reference. And the relevant publications are:

K.N. Chaudhury, D. Sage, and M. Unser, "Fast O(1) bilateral  filtering using
trigonometric range kernels," IEEE Transactions on Image Processing, vol. 20,
no. 11, 2011.
K.N. Chaudhury, "Acceleration of the shiftable O(1) algorithm for bilateral filtering
and non-local means," arXiv:1203.5128v1.

I assume the Gaussian kernel is position shift-invariant. So normalization means the sum of elements inside the kernel equals to 1. For a 3*3 Gaussian kernel, it is

1/16  2/16  1/16
2/16  4/16  2/16
1/16  2/16  1/16

The bilateral filter can also be normalized with the K factor of integral of the product of  Gaussian kernel and intensity kernel. You can try both to observe the effect of normalization on your image data.


A check on the definition of the continuous periodogram, and does it assume ergodicity somewhere? Can someone verify my understanding of what the continuous periodogram is/means, and please tell me if I say something wrong:

As I've learned so far, the power spectral density of a wide-sense stationary stochastic process $x(t)$ is given by

$$S_{xx}(f) = \lim_{T\rightarrow\infty}\frac{1}{T}\mathbf{E}\Big\{|X_T(f)|^2\Big\} \tag{1}$$

where $X_T(f)$ is the Fourier transform of a single truncated sample path. The sample path itself is over $\{{-\infty<t<\infty}\}$, but the truncation is over time $T$. Finally, the expectation is over all sample paths in the ensemble.

1) Is it correct to write then that: if in practice we only have access to a single member of the ensemble which is over a finite observation time $T$, then the simplest possible estimator, $\hat{S}_{xx}(f)$, of the true power spectral density, $S_{xx}(f)$, is found by dropping both the $T\rightarrow\infty$ limit and the expectation (because we don't have any other members of the ensemble), and so

$$\hat{S}_{xx}(f) = \frac{1}{T}|X_T(f)|^2 \tag{2}$$

and that this is known as the continuous periodogram?

2) I have seen some resources which say that the periodogram assumes ergodicity - is this true? If so, where and how exactly does this assumption become necessary? As far as I can see, if we simply drop the expectation in (1) (because we only have access to a single realisation), then there should be no reason to need to invoke ergodicity (which would somehow involve replacing the ensemble expectation with a temporal average).

3) Finally, I would like to ask about the variance. I see that the main problem with the discrete periodogram is that it's variance does not decrease as the number of samples is increased. Is there an equivalent problem with the continuous periodogram, and what is the analogy in that case to be able to calculate the asymptotic variance (since we don't have any discrete points to take $N\rightarrow\infty$). Should it be the time window $T$, and one should be able to show that the variance becomes constant as $T\rightarrow\infty$?

Eq. $(2)$ is indeed the periodogram of the truncated signal where said signal is the sample path $x(t)$ of the random process. There is nothing random about this sample path, which is why dropping the Expectation operator makes sense when one goes from $(1)$ to $(2)$. Dropping the limit is asking everyone to take your word for it that the $T$-second observation interval $\left[-\frac T2, +\frac T2\right]$ is adequate to capture enough of the general behavior of the process. Would increasing $T$ to $T^\prime > T$ make a significant difference in $\hat{S}_{xx}(f)$, at least as far as the central lobe is concerned?  We don't know but we trust you to make an informed judgement.

With regard to your query $2)$ as to whether the periodogram assumes ergodicity, the answer is No and Yes.  No in the sense that Eq. $(2)$ by itself says nothing about the power spectral density $S_{xx}(f)$ all by itself: the right side of $(2)$ is the periodogram of the observed finite-length segment of the sample path of the process. Yes in the sense that it is you who can be said to be assuming ergodicity when you assert that $\hat{S}_{xx}(f)$ is a good estimate of $S_{xx}(f)$.  In general, the averages of the sample paths (e.g. the limit as $T\to \infty$ of the right side of $(2)$) need not have the same properties as the averages over ensembles (e.g. $S_{xx}(f)$), and assuming that they are the same not only means that you are assuming ergodicity (which really deals with limiting or asymptotic values of the averages) but also that $\hat{S}_{xx}(f)$ is pretty close to $\lim_{T\to\infty}\hat{S}_{xx}(f)$, that is, the short observed segment of the sample path is generally representative of the process.


Definition of sampling using delta or indicator function? I just came from a class where the professor showed a slide with the definition of sampling:

But I do not understand how we can multiply a signal $x(t)$ with the delta function $\delta(t)$, as the $\delta(t)$ is infinite at $x=0$. So this multiplication would amount to multiply a real-valued function with $\infty$ which is obviously not what is meant here. I think.

So would it make the formula formally correct if we would swap the $\delta(t)$ with the indicator function?

$$\mathbf{1}_0(x) =
\begin{cases}
1 &\text{if } x = 0, \\
0 &\text{if } x \neq 0.
\end{cases}$$

that equation following "Conceptual scheme" is misleading.  the $\implies$ arrow is not an equal sign.  this is what is true:

$$\begin{align}
x(t) \times \sum\limits_{n=-\infty}^{\infty}\delta(t-nT_s) &= \sum\limits_{n=-\infty}^{\infty}x(t)\times\delta(t-nT_s) \\
&= \sum\limits_{n=-\infty}^{\infty}x(nT_s)\times\delta(t-nT_s) \\
&= \sum\limits_{n=-\infty}^{\infty}x[n]\times\delta(t-nT_s) \\
\end{align}$$

because we simply define $x[n] \triangleq x(nT_s)$.

the difference between the dirac delta and the indicator function is that one integrates to an area of 1 and the other integrates to an area of 0.


Why use both packet detection and symbol timing offset in WLAN 802.11? In IEEE 802.11 Wlan, at the receiver side, we first use packets detection algorithm to detect the start of the packets, and then uses fine symbol timing estimation algorithm to find the start of the data symbol.
However, if we can detect the start of the packets, we should be able to find the start of the data symbol since the length of the preamble is known at the receiver, so why would we use still employ the symbol timing estimation algorithm to find the start of the data symbol. The packet structure:

I guess it because the multipath fading or symbol time offset, but I can't figure out why, could any body explain this? I will really appreciate that, thanks in advance.

However, if we can detect the start of the packets, we should be able to find the start of the data symbol since the length of the preamble is known at the receiver, so why would we use still employ the symbol timing estimation algorithm to find the start of the data symbol.

The start of packet detection isn't accurate enough. When you look into the properties of OFDM symbol timing recovery methods. you'll notice that they use mathematical tricks to achieve higher timing accuracy.
The main function of the packet detection is to know that you need to synchronize, not to synchronize. (and of course, to avoid collisions)


Fast/efficient way to compute Laplacian edge enhancement filter I would like to implement a somewhat smarter Laplacian edge enhancement convolution. Right now it is implemented as (generic 3x3 convolution):

sum = 0
for k = -1 to 1 do
   for j = -1 to 1 do
      sum = sum + h(j +1, k + 1)*f(x - j, y - k)
   end for
end for
g(x, y) = sum

Looking for way to improve to implementation all I could find was separable filter optimisations (as in). However Laplacian edge enhancement is not separable:

octave:1> m = [-1 -1 -1 ; -1 8 -1 ; -1 -1 -1];
octave:2> rank(m)
ans =  2

Is there any litterature I could be missing ? When looking at the coefficients it feels as if somewhat could be done for matrix of type (CENTER=8, EDGE=-1):

E  E  E
E  C  E
E  E  E

OK here I will put an old trick to reduce number of multiplications dramataically, whether it will render a faster result depends on your overall computational architecture of course.

Now as your filter is highly structured, let us denote two values
E and C to describe its elements. Lets decompose your filter into two parts:
$$h[n,m] = h1[n,m] + h2[n,m]$$ such that:
$$h1[n,m] = E*ones(3,3)$$ and $$h2[n,m]=(C-E)*\delta(n,m)$$

here I put a pseudo-language version of your convolution core, that takes advantage of this simple decomposition:

for (n,m)= beginning:end   // it is a double loop
   sum = 0;
   for (j,k)=-1:1         // it is the core double loop
      sum = sum + f(n-j,m,k);   % Here there is no multiplication, saving 9  multiplies
   end:j,k

   g(n,m) = E*sum + (C-E)*f(n,m); % here is two multiplies
end:n,m

Note1: please take care of indices especially at the edges, since I didn't.
Note2: whether your optimizing compiler can already convert your code into this structure, from your previous code, depends on your filter definition. If your filter is defined only by a pointer without its content, the compiler will have no chance to decompose it. But if the filter is statically defined at compile time, then at least in principle the compiler has a chance to analyse its content. Whether they do it or not is another concern...


Quadrature modulation for FM and AM This might be a naïve question, are these correct statements?

Both amplitude and frequency modulated radio signals nowadays use quadrature modulation and demodulation as a mean to transfer and receive radio signals.
Whenever quadrature modulation is used the sampling rate should be equal to bandwidth.

Both amplitude and frequency modulated radio signals nowadays use quadrature modulation and demodulation as a mean to transfer and receive radio signals.

This question is ill-defined – a signal doesn't use quadrature modulation/demodulation, transmitters/receivers do.

So, your question is, if we try to "rescue" it, is

  Do modern FM receivers/transmitters use quadrature modulators

No, not generally. When I receive FM audio with my USRP, I might do that, but usually, I will use a trick and actually do a low-IF sampling.

When a car radio receives signals, it definitely doesn't – simply because IF receivers typically still have technological advantages over direct conversion /IQ receivers, mainly due to lack of DC offset, and because you can get high-Q analog bandpass filters easily. The bandwidths used for FM audio are so laughably small that any modern 16 or 32 bit microcontroller can usually deal with the IF sampling rate.

Now, assuming this is about FM audio broadcasts only: I'd suspect that you'll actually see the different components of the signal first being generated digitally, then digitally mixed to different IFs, and then generated using IQ DACs and mixed up using a Quadrature modulator. But: Broadcasting equipment often has seen quite a few winters, so this might not always be true!

  Whenever quadrature modulation is used the sampling rate should be equal to bandwidth.

No, not generally. Nyquist says that in the complex sampling case, it is sufficient to sample with nothing more than the signal bandwidth, but in fact, you rarely do that, because you then wouldn't have the chance to suppress noise due to oversampling, or to have "spare" bandwidth to do frequency correction. Also, there might be simple technical problems ("can't get an oscillator running at 13.13131313 MHz") that make it impossible to sample at exactly the bandwidth. And,  there's also applications where undersampling is employed.


Smoothing an LFO Sawtooth Wave I'm trying to smooth a Sawtooth wave coming out of my LFO.   Effectively the transition from 1.0 to -1.0 in a 1 sample step is causing clicks and pops if I modulate my audio signals amplitude with this LFO output.

What is the best way to smooth this basic waveform to eliminate the clicks without adversely affecting performance (significantly)?

Reason that it's clicking is that you are jumping from 1 to 0 in the time that the audio signal reads the last of the 1 control signal and then the 0 control signal.

You have to slow fast signal transitions over multiple samples, for example 1 to 0 should take 5ms or 10ms which is dozens or 100ds of samples.

To do that, you can use a lowpass filter at about 300Hz more or less, which can actually be very cheap for a single monophone LFO its less that 1% cpu useage, or you can use a signal smoothing algorithm, for exampel a linear interpolation function or a SmoothStep which is awesome general knowledge for audio programmers.
http://paulbourke.net/miscellaneous/interpolation/

BTW what is the platform? i can't give you code exampels if i dont know what code you use.

the cool thing about smoothstep is how easy it is, its very good knowldedge in synths:
https://en.wikipedia.org/wiki/Smoothstep

https://documentation.apple.com/en/logicstudio/instruments/index.html#chapter=A%26section=3%26tasks=true

check the filtered sawtooth there.

You can also a smoother that says "if difference between last signal and this signal was larger than 0.02, signal - 0.02... and create a secondary smothed signal... there are fairly simple smoothing algo's based on that it's pretty straightforward.

There is also an entire page on the topic.


Mathematical operations with given discrete frequency response function Let's say I have a given discrete frequency response function H(w) and corresponding n frequencies as MATLAB arrays, ranging from 0 to 512 Hz.

If I do an n-point FFT on a discrete time domain signal x(k), how do I filter this frequency domain signal with given FRF? Is it just element-wise multiplication of the FFT array with the FRF array, followed by an IFFT to get the filtered signal y(k) in time domain or is it more complicated than that?

To take it one step further, what if I only had a filtered signal y(k) and my discrete FRF, could I compute the original signal x(k)?

Is it just element-wise multiplication of the FFT array with the FRF
  array, followed by an IFFT

No. Filtering in the frequency domain is not trivial since it's equivalent to circular convolution (not linear) and will result in time domain aliasing. Google "overlap-add" or "overlap save" for a description of adequate algorithms

  To take it one step further, what if I only had a filtered signal y(k) and my discrete FRF, could I compute the original signal x(k)?

Sometimes. In essence what you are trying to do is an "inverse filter". This is only possible if the filter itself is invertible. In other words: information that's been lost in the filtering process cannot be recovered.  In practice that means you can't have any zero (or too small) amplitude in your filter and you need decent signal to noise at all relevant frequencies


How to Detect a Inhomogeneity Region in Image I have an image that includes object and background. However, the object appears some inhomogeneity region due to illumination. My work is that how to detect inhomogeneity region. Which is feature can represent it? Example: low inhomogeneity-high/low feature value....Thank you so much

or

There are many properties of inhomogeneity:

Local Variance / STD.
Local Histogram.
The Gradient Function
Histogram of the Gradient.
Mean versus the Median / Mode.


issue in DFT of OFDM by blocks and N-by-N I'm trying to study and understand the OFDM, and during coding the OFDM in MATLAB, I faced this issue:

suppose we have the Numsymbol = 256, and dataqpsk our used data, The issue which I face that when I perform iFFT block by block, that seems to be ok, but when trying to make it N-by-N .. that gives a problem.

here is the steps. now I will make the iFFT by blocks :

clear all; clc;

N_carrier = 256;  % number of carrier
N_symbol = 256;  %number of symbols
N_FFT = 256;  %length of FFT
N = N_carrier;
L_cp = N/4;   % the CP

dataqpsk = rand(N_carrier,N_symbol);  %Example for modulated data
size(dataqpsk);
for isij=1: N_symbol
    % input modulated data
    dataqpsk1=dataqpsk(:,isij);
    size(dataqpsk1)
    D_map = dataqpsk1;

    %% doing such mapping
    D_IFFT_fore=zeros(N_FFT,1);
    D_IFFT_fore(1:(N/2)) = [D_map(((N/2)+1):N)];
    D_IFFT_fore((N_FFT-((N/2)-1)):N_FFT) = [D_map(1:(N/2))];
    D_IFFT_aft=N_FFT.*ifft(D_IFFT_fore,N_FFT);%×÷ifft

    %add GI and interpolation.
    D_cp=[D_IFFT_aft(N_FFT-L_cp+1:N_FFT);D_IFFT_aft];
end

But now when I try to make it as DFT N-by_N I don't get, as below:

clear all; clc;

N_carrier = 256;  % number of carrier
N_symbol = 256;  %number of symbols
N_FFT = 256;
N = N_carrier;
L_cp = N/4;  %CP

dataqpsk = rand(N_carrier,N_symbol);
size(dataqpsk);
D_map = dataqpsk;
D_IFFT_fore=zeros(N_FFT,N_FFT);
D_IFFT_fore(1:(N/2)) = [D_map(((N/2)+1):N)];
D_IFFT_fore((N_FFT-((N/2)-1)):N_FFT) = [D_map(1:(N/2))];
D_IFFT_aft = N_FFT.*ifft(D_IFFT_fore,N_FFT); % ifft
D_cp(1:L_cp,1:N_FFT)=D_IFFT_aft(N_FFT-L_cp+1:N_FFT,1:N_FFT);

The problem is when trying get the results of DFT N-by_N, I get from the step of D_IFFT_fore(1:(N/2)) = [D_map(((N/2)+1):N)];  a matrix of Zeros !!

please, any help?

thank you

Thank you Dear msm, and here is the complete code:

clear all; clc;

N_carrier = 256;  % number of carrier
N_symbol = 256;  %number of symbols
N_FFT = 256;
N = N_carrier;
L_cp = N/4;  %CP

dataqpsk = rand(N_carrier,N_symbol);
size(dataqpsk);
D_map = dataqpsk;
D_IFFT_fore=zeros(N_FFT,N_FFT);
D_IFFT_fore(1:(N/2),:) = [D_map(((N/2)+1):N,:)];
D_IFFT_fore((N_FFT-((N/2)-1)):N_FFT,:) = [D_map(1:(N/2),:)];
D_IFFT_aft = N_FFT.*ifft(D_IFFT_fore,N_FFT); % ifft
D_cp(1:L_cp,1:N_FFT)=D_IFFT_aft(N_FFT-L_cp+1:N_FFT,1:N_FFT);

thanks again

Try using something like the following. You want to map the rows this way for all columns, hence the colon to denote all columns. Double-check that you get what you want from these. I didn't run it to verify it gives exactly what you want, but it should fully populate all columns, I believe.

D_IFFT_fore(1:(N/2),:) = D_map(((N/2)+1):N,:);
D_IFFT_fore((N_FFT-((N/2)-1)):N_FFT,:) = D_map(1:(N/2),:);


How is the number of amplitude levels decided in modulation techniques? For M-QAM modulation, where e.g. $M = 64$ there are 8 amplitude levels. Each level is represented by a symbol that is further represented by bits through the relation:
$$b=\log_2(M)$$
This means that each symbol is encoded by $b$ bits.

I don't quite understand where 64 comes. Does it mean that the imaginary and real component are included ($8\cdot 8$) ?
For other modulation methods such as 16-QAM, where is the 16 coming from?

There are $M$ symbols in an $M$-QAM constellation where two carriers I and Q are amplitude modulated, and there are $N=\sqrt M$ amplitudes for each I and Q carriers (for $M\ge4$). In $16$-QAM there are 16 symbols on a $4\times4$ constellation. That is, there are $4$ amplitude levels for the I and $4$ levels for the Q component. Of course, $\log_2M$ bits are conveyed by each symbol. Since $M$-QAM normally uses carriers with symmetric amplitudes, $M$ is usually a power of $4$. For instance, $M$ can be $4, 16, 64, 256,\cdots$.


Analyzing a 16-QAM modulated waveform? Okay I'm really confused about analyzing the  waveform I generated using 16-QAM on MATLAB... I transmitted a $y = [0 \ 1 \ 2 \ 3 \ 4 \ldots \ 15]$ discrete signal (in order) to qammod function in MATLAB usign 16-QAM. Its constellation diagram is the standard 16 symbols on the combination of values of $-3,\ -1,\ 1,\ 3$ on Q and I. Then I modulated the I and Q components respectively with a cosine and sine, with a carrier frequency of $36\textrm{ kHz}$. My result is this:
I know it follows a complex envelope, but how do I know that I got a correct waveform based on my original signal?

Look at the magnitude and angle of the signal relative to your carrier and you should see the same thing you would see if you looked at your angle at baseband using the I and Q samples.

It is difficult to see anything with the carrier there (and I would never simulate with an actual carrier when I can do everything at baseband, just in case you put it at the carrier for simulation purposes). To inspect this waveform, I suggest rotating it to baseband by multiplying it by $e^{-j2\pi f_c t}$, low pass filter the high frequency component assuming your input was real, and then plot the complex signal on an I Q complex plot. (Real on one axis and Imaginary on the other). You should expect to see the signal moving through your 16QAM constellation points with a trajectory that depended on the pulse shaping filter used. Alternatively you could also look at the real and imaginary separately in an eye diagram.

Here are some plots I have showing that may help illustrate the different cases:

On the left is an eye diagram and the associated constellation diagram for QPSK, with two cases for the IQ (real-imaginary) constellation; one where it was derived from an incorrect sampling point (timing error) and the other where each symbol was properly sampled, as indicated by the vertical bars in the eye diagram. Note the eye diagram is simply our waveform versus time over one symbol period, plotted again and again (as in infinite hold on a scope) such that we see every possible trajectory from one symbol to the next, which is dependent on the history of all prior symbols over the memory of our channel.

On the right is the constellation diagram for 16QAM, where the large one is what it would look like with a carrier offset (similar to your case), such that it is essentially our 16 points spinning (which doesn't tell us much!), and the other is after removing the spin with $exp^{+j2\pi f_c t}$ or $exp^{-j2\pi f_c t}$ depending on which way the diagram is rotating. In your case with a real signal both spins are present, so multiply by either and low pass filter the resulting higher spin portion. One direction will have an inverted spectrum compared to other, so if you go the wrong way you may have an inverted result, but this is easy to undo by simply swapping I and Q).  The smaller plot on the right is the resulting 16QAM constellation plotted as real vs imaginary of ONLY the correct sampling locations within each symbol were chosen (no timing error as in the QPSK case shown), but we do see evidence of phase error and some amplitude error (hence the spreading of the dots), as well as the possibility of other distortions such as multipath etc causing intersymbol interference; all having the effect of increasing the size of the dots in our constellation (leading to errors!).

Below is another example for QPSK showing the constellation and eye diagram side by side, and also showing the significance of sampling at the correct sample positions (in red) versus the overall trajectory of the waveform.

And finally here is an example eye diagram for 16 QAM:


Standard deviation of Gaussian corresponding to perfect low pass filter [1,1] For a time discrete signal, the rect-filter (for 1D signal [1,1]) is a perfect low pass filter. I was wondering, what is the best gaussian approximation of this filter?

For a current problem, I cannot use binomial filters but have to use Gaussian filters for smoothing. Ths question is thus, what is the standard deviation of a Gaussian filter best approximating the [1,1] rect-filter. 'Best' could e.g. mean is has similar probperties in the frequency domain.

Suppose you want a low pass filter supported on $f \in [-1,1]$. The Fourier transform of a Gaussian is a Gaussian. So let's assume our Gaussian filter has a frequency response of the form $G_\sigma(f) = \frac{1}{\sigma \sqrt{2 \pi}} e^{-f^2/\sigma^2}$. For the best Gaussian fit in least squares sense, let us minimize the squared error ($\mathsf{SE}$):
$$
\mathsf{SE}(\sigma) = \int_{-1}^1 |G_\sigma(f)-1|^2 \;\;\; d\,f
$$
by choice of $\sigma$. This calls for some numerical evaluation.

I used Mathematica:

Integrate[ (1/s/Sqrt[2 Pi] Exp[-f^2/s^2] - 1)^2, {f, -1, 1}]

which gave me

2 - Sqrt[2] Erf[1/s] + 1/(2 Sqrt[2 Pi] s) Erf[ Sqrt[2]/s ]

Then plot the result as a function of $\sigma$ (i.e. s).

I see a minimum around $\sigma =0.7$.


follow up to previous question " is this an impulse response or step response" This question is a follow up to the gorgeous explanation of how the leaky integrator model becomes  a differential equation. The link to that question and answer is below.

should this be viewed as an impulse response or step response

I have two questions related to the explanation. Hopefully Fat32 will not mind answering but anyone else is great also.

1) At the bottom of the first page, the author has  2 statements that start with "such that".

The beginning of the first one is $\delta(t) = x(t)  \cdots $  and the beginning of the second one is  $u(t) = x_(t) \cdots $

Clearly, the equality sign in those statements does not mean the usual equality. Could  that notation be clarified because I don't follow it ?

2) At the very end of the answer, the author derives the continuous time differential equation for the leaky integrator model:

$y^{\prime}(t) + \alpha \times y(t) = x(t)$.

In the link I was originally looking at, the differential equation is discrete rather than continuous but aside from  that, the same equation as that immediately above. ( proportionality constant is $1/k$ rather than $-\alpha$ ). Of course, the discrete version of above would be

$ dy_{t} = (- \alpha \times y_t + x_t ) dt $

Where the subscripting denotes the discrete unit of time for $t = 1, \ldots n$ and $dt$ denotes the width between these discrete units.  This discrete version above is what Lesson 12 uses to generate the plots for how the level of the water in the bucket changes.

So my question is the following:  Would this discrete equation directly above  just be derived by discretization of what was derived by Fat32 or would the leaky integator model derivation need to change somewhere to justify discreteness ?  Clearly, the notion-concept of $exponential declining memory goes away in the discrete version of the model so that makes me think that maybe a whole new derivation is required for justification ? If it is required, then I imagine it is somewhere but I can't find it. Thanks a lot.

                                                Mark

Let me provide an answer,

1-) In the answer I have used the notation $x(t)=\delta(t)$ and $x(t)=u(t)$ to indicate that the input to the system is a dirac impulse and unit-step function respectively. The correct notation should be $x(t) \equiv  \delta(t)$ which means identical, by definition. Triple equality sign is almost always omitted and the double equality sign is used instead in signal processing context. (it's therefore not asking for the solution of an equation, which is generally understood from double equality sign).

2-) The subject of differential and difference equations is rich. According to the usual literature, continuous time-systems are described by differential equations and discrete-time systems are described by difference equations. When a continuous-time signal is sampled and the associated continuous system is transformed into a discrete-time system the corresponding differental equation is also discretized and converted into an equivalent difference equation. One application of such discretization is the usual Kalman filering equations.

I've analysed the continuous-time leaky integrator model whose input-output relationship was assumed by definition as $y(t) = \mathcal{T} \{x(t) \} = \int_{-\infty}^{t} e^{-\alpha(t-\tau)} x(\tau) d\tau $. Given this definition of the system the corresponding continuous-time LCCODE is obtained as in the answer. In the web-link only the differential equation relation is provided as $ dy/dt = s - y(t)/k $ where $y(t)$ is the water level (the output) $s=s(t)$ is the input and $1/k$ is the leakage constant, the differential equation is then $$ \frac{dy}{dt} + \frac{1}{k} y = s $$ or equivalently $$ \frac{dy(t)}{dt} + \frac{1}{k} y(t) = s(t) $$ This differential equation is the same as $ y' + \alpha y = x$ after replacing $\alpha = 1/k$ and $s(t) \equiv x(t)$.

On the other hand in your original question, you refer to a discrete differential equation that I cannot clearly understant whether it's a difference equation or a differential equation (the terms discrete and differential used together as a mix) or it's just another application of discretization as I defined above.

However the following excerpt from the original web-reference indicates to me that:

  We can easily simulate this leaky integrator in discrete steps of time
  by changing the value of 'y' on each step according to the equation
  above:

There is only a continuous-time differential equation and it's solution is numerically computed by the software, which is indicated as the simulation of the ODE using discrete steppings (the numerical computation)

So as I understood given the continuous time differential equation
$$ y'(t) + \alpha y(t) = x(t) $$ one can reach the numerical version as
$$ \frac{d y}{dt} = -\alpha y(t) + x(t) $$
$$ dy = ( -\alpha y_t + x_t) dt $$

Where $y(t)$ is replaced by $y_t$ just as a matter of notation. So they are the same thing. The web-link numerically solves the contiuous differential equation for various input stimulus.


Relationship between SNR and link speed I've been trying to model some radio signal properties of wireless networks wrt its usefulness. For example, most WiFi network adapters would decrease their transmission rate if the overall quality of the signal is bad so as to avoid re-transmission, errors, etc.

What is the relationship of the received signal strength (measured in dBm) and the signal-to-noise ratio (SNR) wrt the overall signal quality? Ideally, I am looking for an equation or paper reference.

Thanks!

Signal strength is only important in determining signal quality insofar as it affects the signal-to-noise ratio (SNR).  As the name implies, SNR is determined by dividing the signal power by the noise power in the same bandwidth ($\frac SN$).

Claude Shannon, the pioneer of information theory, came up with a formula that gives an upper limit to how much information (you can think of this as bits) can be transmitted for a given SNR and bandwidth $B$ (I will note, though, that Wikipedia calls it the Shannon-Hartley theorem).  The formula is as follows-

$$
C = B\log_2\left(1 + \frac{S}{N}\right)
$$
where $C$ is the channel capacity in bits/s, $B$ is the bandwidth, and $\frac{S}{N}$ is the SNR.  Modern transmitters/receivers try to get as close to this limit as they can through a number of means.  One way is to increase the bits per symbol as the SNR improves (e.g. going from QPSK to 16-QAM).  Another is to reduce the error correction coding as the SNR improves, and increase it when it worsens.


does scrambling increse the bandwidth In LTE systems the in the scrambling section, a gold sequence is used. But I've studied that gold sequences are used in spread spectrum communication. This generally increases the bandwidth. So doesn't this mean that scrambling results in an increase of bandwidth.?

regards,
phani tej

No, scrambling does not increase the bandwidth.  Scrambling does even out the bandwidth, such that the power is spread evenly over the signal bandwidth.

You are thinking of chipping, which is used in DSSS systems.  Scrambling and chipping are different processes.


1D first-order derivative mask I'm currently reading about segmentation, and have come across finding edges in an image. I then have this image:

Now, I understand how the mask used for the second-derivative works, and give the values for that, but I can't figure out how the mask for first-order derivative works ?
The second-order 1D mask is just [1 -2 1], but how does the first-order one look like ? If there even is anyone. Of course I do know, that the first-order derivative is just the slope of the line, and from that it's easy to see what the values should be. But mask wise, I'm not sure how it works.

Thanks in advance.

As you can see in the following image, the image shows the first order 1D derivative.
Now you can write this equation in terms of the previous pixel rather than the following pixel.

For a 1D differentiation, you are only interested in either the x direction (Horizontal) changes of pixel intensity values or the y direction (Vertical).

Therefore the mask will be [0 -1 1] in case of next intensity value.
Or [1 -1 0] for previous intensity value.

You can apply the same idea for vertical differentiation by taking the transpose.


Noise when dividing or multiplying signal in C I'm trying to create a echo effect and have managed to add multiple delays but when i try to divide or multiply these numbers the signal is completely distorted and adds a lot of noise to the system.

This is my some of my program written in C.

            unsigned char txBuf0[AUDIO_BUF_SIZE];
            unsigned char rxBuf0[AUDIO_BUF_SIZE];

            int i, n, x;
            unsigned short Stor0[2000], Stor1[2000],Stor2[2000], Stor3[2000], z, ze;

            for(i=4; i < AUDIO_BUF_SIZE; i+= 4)
            {
                x= i/4;
                if(lastFullRxBuf == 0){
                    Stor3[x] = Stor2[x]/2;
                    Stor2[x] = Stor1[x]/2;
                    Stor1[x]= Stor0[x]/2;
                    Stor0[x] = 0;
                    Stor0[x] =((rxBuf0[i-4]) | (rxBuf0[i-3])<<8);

                    z = (Stor0[x]) + (Stor1[x])+ (Stor2[x]) +(Stor3[x]);

                    rxBuf0[i-4]= z;
                    rxBuf0[i-3]= z>>8;

                    txBuf0[i-4]=rxBuf0[i-4];
                    txBuf0[i-3]=rxBuf0[i-3];
                }

In processing integer audio samples, simple non-scaled integer division can introduce quantization noise (rounding errors).

You might try doing your amplitude scaling using a higher resolution format (floats or scaled integers with more fractional bits) for your arithmetic, and then perform appropriate noise filtering (dithering or noise shaping of the remainder from the quantization process) if a shorter integer result is required.

Also watch out for integer overflow (e.g. make sure your number formats can hold the full range of expected arithmetic results).


Energy Preserving Laplacian Pyramid Both Discrete Wavelet Transform (DWT) and Undecimated DWT possess an important property of energy preservation: on each level

$$\sum_i W_i^2 + \sum_j V_j^2 = \sum_k X_k^2$$

where $W$ and $V$ are detail and approximation coefficients, respectively.

It turns out that Laplacian Pyramid (LP) generally lacks this (important) property. Is it possible to construct a filter for LP decomposition that would preserve energy? I'm ready to sacrifice the simplicity of backward transform.

The bad news is that I disagree with the first assertion: there exist non-orthogonal DWT (biorthogonal for instance).

But the good news is that this is not important here, and the answer is yes. A Laplacian pyramid is an overcomplete expansion. It can not be treated as a basis, but as a frame. And the concept akin to orthogonal bases in frames is the tight frame.

You can get details on such constructions in Framing Pyramids (2003), whose abstract is:

  In 1983, Burt and Adelson introduced the Laplacian pyramid (LP) as a
  multiresolution representation for images. We study the LP using the
  frame theory, and this reveals that the usual reconstruction is
  suboptimal. We show that the LP with orthogonal filters is a tight
  frame, and thus, the optimal linear reconstruction using the dual
  frame operator has a simple structure that is symmetric with the
  forward transform. In more general cases, we propose an efficient
  filterbank (FB) for the reconstruction of the LP using projection that
  leads to a proved improvement over the usual method in the presence of
  noise. Setting up the LP as an oversampled FB, we offer a complete
  parameterization of all synthesis FBs that provide perfect
  reconstruction for the LP. Finally, we consider the situation where
  the LP scheme is iterated and derive the continuous-domain frames
  associated with the LP.

And yes, the pseudo-inverse can be less simple, but you will gain more optimality in the presence of noise.

Notably, such design allowed the construction of contourlets, a multiresolution directional tight frame designed to efficiently approximate images made of smooth regions separated by smooth boundaries.


Amplitude and Phase spectrum of a signal Please help me to solve this question.

When the spectrum is complex conjugate symmetric around $f_s/2$ (and multiples etc.) the corresponding signal is real valued and vice versa.

This holds for your example, meaning your signal should be real-valued. Your phase spectrum is interesting, since zero phase shift in all frequencies means your spectrum is real valued as well.


How to get the X and Y Coordinates of the circles from an image? I am building an simple website to show the first train arrival and last train departure timing in Singapore. I already wrote a simple code which could be found here. I am placing an circle on top of each station so when hovered
it will display the timing for each station.

Thus far I found the location(X and Y Coordinates) of each station manually added it to the code. But then I just realized that the map I am using is not the latest MRT Map. I don't want to go through the process again waste time for the new map. I am looking for a easy way to do this from the code.

There are two techniques I can think of right now:

Find Equations for each Line assuming circle as points.Then place
the circle at the x and y coordinate. I understand  that one only
Harbourfront to Punngol(Purple) line is the easiest because its a
straight line. I am not so sure what is the equations for other
lines.
Use template matching technique. Take a copy of the circle and look
for that circle in the whole image and get the coordinates.I
understand OpenCV can do template matching but how to do it in
JavaScript? Is there any simple JavaScript  image library to do
template matching?

Any other suggestions?

It is indeed possible to detect the center of the circles by various image processing techniques.

For example, convert to grayscale, apply dilation until the edges disappear, then connected components and compute the centroids.

But this only gives you points, in no specific order, and this might be of little help to you.

It is possible to match the new map against the old one, match the circles, to transfer existing info, but you'll have to check the correct transfer "by hand", anyway.

I wonder why this should be done in JavaScript, as I think that this is an offline operation.


How to add repeated beeps to sound created with a pause in between each beep in MATLAB I have 2 sounds that I created and trying to make them beeps. Maybe 4 or 5 repeated beeps but have a pause in between each beep but make it still part of the signal. I think I have it so that the sound is a single frequency but just need make it a repeated with pauses and still part of the same signal. Any help would be appreciated I have the code listed below.

Fs = 44100;
N = 400;
t = (1:N)*(1/Fs);
freq = 600;

sound_array = zeros(N, 2);
%starts N x 2 matrix
sound_array(:,1) = sin(2*pi*freq*t);
sound_array(:,2) = sin(2*pi*freq*t);

%Left and Right channels

sound_L = sound_array;
sound_L(:,2) = 0;
sound_R = sound_array;
sound_R(:,1) = 0;
plot(sound_array);

%Sound 1
sound(sound_L, Fs);
%Sound 2
sound(sound_R, Fs);
%Play consecutive
sound(sound_array, Fs);

Any pause can be done by using vector/matrix concatenation. Let's say you want to make a signal with pause of length P after it, then do the following:

P = 400;
beep = [sound_array; zeros(P, 2)];

Now if you want to make arbitrary number K of your beeps then use the repmat function as follows:

K = 5;
beeps = repmat(beep, K, 1);

Please read more about matrix concatenation here if you wish.

I believe that with your current code, because you use sound function, signals will be played at the same time. I suggest you to use the pause function to separate playbacks, i.e.:

sound(...)
pause(1)
sound(...)

To create a 1 second pause between playback starts.


Why exactly is the Guassian filter a low-pass filter? In image processing, a Guassian filter is used to blur an image. A low-pass filter is a filter that attenuates the high frequencies, preserving only smooth variations in the provided image. Is a Guassian filter a low-pass filter because it smooths the edges (or, in general, the parts of the image where there is a high variation of the intensity)? Why exactly is the Guassian filter a low-pass filter?

Because the Fourier transform of Gaussian in the spatial/temporal domain is also a Gaussian in the frequency domain. Since a Gaussian drops off rapidly with frequency, it's a low pass filter.


Phase measurement with FFT for DMM Keithley 7510 I am currenly working on a measurement driver for the DMM Keithley 7510. I implemented a harmonics measurement using FFT (from MathNet library).

First I simply sample my input voltage signal and apply the FFT.
From the complex values I get the amplitudes and the phases of the harmonics. The phase results are the phases relative to the fundamental.

I verified my implementations for the amplitude and phase measurement functions by feeding the functions with an "ideal" signal - a sine signal with harmonics generated in MATLAB.

The amplitude measurement works perfectly fine, but I have huge problems with the phase.

Firstly, with higher frequencies, the deviation of the phase increases. This is due to the input characteristics of the DMM, which is not the main problem. I did my measurements mainy at 50Hz.
Secondly, the phase depends on the number of samples. I did not analyze this too much, because I have a fixed sample rate for my measurement driver.
Now the main problem: there seems to be a correlation between the phase and the amplitudes of the harmonics (including the fundamental). This means that by changing the amplitude of a harmonic, the phase results change for all the harmonics. For instance an increase does not necessarily mean an increase in phase: for some harmonics, the phase increases while for others it decreases. I could not find any reasonable explanation for this behaviour.

I tried many different ways to solve the phase problem:

Windowing: I thought, that it might somehow help. Unfortunately, I found pretty much no sources on the effect of windows on the phase measurement.
Power of 2 for FFT: I have a power of 2 number of samples on which I apply the FFT.
DFT: I also tried to apply DFT and calculate the complex value at exactly the frequecies where I wanted. There is as good as no difference between FFT results and DFT results.
Amplitude interpolation: I found a paper where I found I formula for the amplitude interpolation. By finding the peak amplitude in the frequency domain, I calculated the corresponding phase, but no success.
Here I later found, that at the peak amplitude of a harmonic in the frequency domain, the corresponding phase is not the expected/ideal phase. This is also something, that I could not understand ...

My question now is whether I forgot any other things I should check for.
Are there any parameters that I should be careful about when computing the phase?

Any suggestions what I could do next?

First off, you should realize that the FFT is just an efficient implementation of the DFT.  Power of two sizes are the easiest for it, but modern implementations have slick tricks to do other sizes.  Leave that up to the library.  The results will be the same.

If your peaks are sufficiently well spaced (two or three bins apart) the following article will show you how to efficiently calculate the phase and amplitude when a peak falls between bins.

Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

For a single pure tone, the results are exact providing you know the exact frequency.  You can find several other articles of mine that will show you how to get a very good frequency estimate, exact in the single noiseless pure tone case.

You will have to read the referenced articles (also mine) to understand how it really works.

Ced

If your signal is a repeating waveform, the best approach is to define your frame on a whole multiple of cycles.  When you do this, each harmonic will fall in its own bin and there will be spacing (depending on the waveform count) from the nearest bins.  When a signal is not a whole number of cycles, that's when the peak falls inbetween the bins and the values get distributed in the nearby bins.  This is called leakage and it causes interference.  The equations I referenced do a pretty good job of ignoring this interference, but to get as good of a reading as you would by framing a whole number of cycles, you need to remove the effects of nearby estimates (subtract them out of the DFT using my bin value formulas) before measuring a peak.

Also, you should be aware that if you are measuring something like a piano, not all the overtones are harmonic.  This prevents you from centering all the bins simulaneously.  See:

https://www.dsprelated.com/thread/7902/the-spectral-complexity-of-a-single-musical-note

The easiest way to really understand a DFT is to realize that each bin index (zero based, not MATLAB one based) corresponds to a tone which has the same number of cycles per frame as the the index number.  Thus if you frame five cycles of your periodic signal, ideally, you expect values in bins 0, 5, 10, 15, etc.  and all the inbetween bins will be zero.  This is independent of the sample count within the frame.  What a higher sample count buys you is more bins at higher frequencies,i,e, you are moving the Nyquist limit up.


Capacity and BER as performance of transmission system We study the performance of the transmission system by channel capacity or BER/SER. If the system performance is studied as BER/SER, does it make sense to study a channel capacity?

Channel Capacity and Bit/Symbol Error rate are two different key performance indicators (KPI) of a system.
BER/SER, as it's obvious from its name, is  the number of Bit/Symbol errors divided by the total number of transferred Bits/Symbols. and you can use channel coding to improve it so it's not a fixed KPI.
while Channel Capacity is maximum rate that the system can achieve and it's a fixed KPI.
so it make sense to study the channel capacity of your system. by studying channel capacity you find maximum rate that can be achieved and by improving BER/SER you can improve rate of system to get close to channel capacity.


ADC sampling rate and signal bandwidth for FHSS If i am correct; for a radio receiver the ADC sample rate will have to be higher than at least twice the BW of the data (disregarding sparse sampling techniques). So for frequency hopping spread spectrum (FHSS),when the receiver is unaware of the transmitter's hopping pattern, does that mean the ADC sample rate needs to be higher than twice the entire hopping bandwidth?
Let's say I use amplitude shift keying, data BW is 500kbps, and the entire hopping bandwidth is 24MHz. If I use a 2Msps (4x data BW) ADC it looks like I will still recover my data; sure it under-samples some of MHz modulated frequency but it does capture entirely my data
correct? if not do I miss something?

So for frequency hopping spread spectrum (FHSS),when the receiver is unaware of the transmitter's hopping pattern, does that mean the ADC sample rate needs to be higher than twice the entire hopping bandwidth?

Yes.

If I use a 2Msps (4x data BW) ADC it looks like I will still recover my data; sure it under-samples some of MHz modulated frequency but it does capture entirely my data.

No.  See above.  It'll alias whatever is in the entire 24MHz band into the 1MHz that it can resolve, or you'll precede it with an anti-alias filter and lose $\frac{23}{24}$ of your signal.
On the bright side, if you're taking 24MHz worth of bandwidth and only using 1MHz of it, that's close to $14\mathrm{dB}$ of processing gain, so while you need an ADC that's 24 times faster, it can be 4 bits narrower.


Feature normalization (scaling) for Hyperspectral images have been straggling with the concept of feature normalization for the hyperspectral images. Respect to my problem I have attached a picture of my problem which clearly stands the issue which I have.

There are different methods for the normalization. If it is for visual graphs, you can divide it into bands which can be defined in realtime, so you'd write a program which presents visual graphs of frequency A->B where A and B can be changed by the user. That helps you see display it graphically. using X and Y axis of the mouse to define A->B bandwidth. Then you can scan backwards and forwards through all the frequencies interactively while also normalizing them, using an amplitude meter beside the graph, to say the overall max amplitude of that bandwidth.

If the normalizing process is for high fidelity, balanced post processing and statistics, then it's best to find the highest value at any frequency and divide all the pixel amplitudes by that high value so that it equals 1.


Eigen Value and Eigen Vectors in Image Processing I have been reading about eigen values and eigen vectors but i haven't be able to find any decent explanation relating their application in image processing/computer Vision for example say face detection and eye detection. Can somebody please through some light on this?

I will try giving you some intuition.
The SVD says each matrix can be decomposed into 3 operations - Rotation, Stretching (Scaling) and the another Rotation.

What matters is which directions are scaled and how.
Directions are vectors (Pointing some direction).

The SVD has many uses in Linear Algebra.
One its most known use is Low Rank Approximation of a matrix.
Low rank meaning to try represent the same data with lower number of degrees of freedom.

It is immediately rings a bell for compression and dimensionality reduction (Also the SVD is practically the basis for Karhunen Loève Compression).

Image Compression

Think of decomposing the image into patches of $ n \times n $ (Where usually $ n = 8 $ or $ n = 16 $).
Now build a matrix of those patches (Each patch as a column).
Then, using the SVD you will be able to find a low rank representation of this image which is basically a compression.
Why isn't it used?
Because this is the most effective compression in $ {L}_{2} $ sense yet the base it uses is adaptive to the image. Namely with each compressed image you'll have also to deliver the base which hurts the practical compression ratio.
But this is the idea.

Eigen Faces

This is a nice implementation of idea using SVD.
You can read about in the great article on Wikipedia - Eigen Face (Including Code).
Basically, again, you can use the low rank representation either to compress the data base of faces and / or use it as a method to recognize faces.

Compressed Sensing

The SVD is heavily used in the Compressed Sensing frame work.
For instance have a look on the Dictionary Learning method called - K-SVD.

The Wikipedia article about the SVD is really great place start reading about it - Singular Value Decomposition (SVD).


Transform Function with Non Linearity I'm a newbie to Signal Processing - my apologies if this question is too obvious (I'm a financial trader trying to use DSP techniques).

For a linear filter:

$y[n] = (1-p) x[n]+p y[n-1]$

we can the write Transfer Function as

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

But instead, suppose I have a non-linear filter:

$y[n] = y[n-1]+(1-p)\frac{1}{2}(x[n]e^{-2 y[n-1]}-1)$

How can I write a Transfer Function ?

For example, if I simply replace $y[n-1]$ by $y[n]z^{-1}$ and rewrite I get:

$H(z) = \frac{(1-p)\frac{1}{2}(x[n]e^{-2 z^{-1} y[n]}-1)}{1-z^{-1}}$

which is a function of $x[n]$ and $y[n]$ in addition to being a function of $z$?

Can this be written merely as a function of $z$?

Transfer function tries to describe a system using the Fourier Transform (More generally the Laplace Transform / Z Transform).
Yet in order for the Transform function to dependent on the frequency and amplitude only the input system must be LTI system.
Hence there is no Transfer Function in the Fourier Domain to your system which can be described the way you want.
Linear Time Invariant (LTI) System
By Convolution Theorem a convolution of LTI System can be described by a transfer function in Fourier Domain.
What do mean by a transfer function? We mean having $ H \left( \omega \right) $ which is independent of the input.
As for any output and input we can calculate the Fourier Transform and have the ratio of them as the Transfer Function. But only for LTI the result of this ratio will be constant.