Carrier Tracking Loop: Tracking noise vs Phase Noise Trade-off I would like to compare two different implementations of carrier tracking loop which differ in the way phase detection happens in the PLL. I am using PSK modulation for the data.
In the first implementation, a mixer is used as a phase detector and the implementation looks like this:

Here, to ensure that the PLL does not track the PSK modulated data, the bandwidth needs to be sufficiently below the symbol rate to ensure that the tracking noise is not too high. This I can understand because if only data (without any noise) is passed through the PLL, low frequency components of data will be tracked and will contribute to noise. This can be seen in the figure below, where the output data from the PLL is high pass filtered because of the Loop:

Another implementation where the phase detector is implemented as a vector multiplier is explained here. In this implementation, I do not clearly see if there would be such kind of tracking noise vs phase noise trade-off. Because if I pass data without noise, than phase detector output will be zero and there would be no filtering of the data by the PLL. The plot for the PLL output with vector multiplier phase detector is shown below and I see that the tracking noise is non-existent even though I chose bandwidth very close to the data rate.
Why would such a trade-off exist for this kind of carrier tracking implementation?

EDIT
For the mixer implementation I am using ideal phase detector:

delta_phi = Kd*(phi_ref - phi_vco) %% Kd is the PD gain

I use complex exponential to implement vector multiplier:

phi_rotation = ref.*vco';  %%% Both and vco are complex exponentials
symbol_estimate = conj(pskDemodulator(phi_rotation));
delta_phi = angle(phi_rotation.*symbol_estimate);  %%% imag would give the same result.

Block diagram Approach 2

Block Diagram Approach 1

The question is specific to optimizing loop bandwidth for the decision directed carrier tracking loop with one sample per symbol. In other posts such as this one PLL for Phase Demodulation and Carrier Tracking I have detailed the consideration for not making the loop bandwidth too low due to increasing contributions from LO phase noise as well as less ability to track system dynamics (rate of Doppler change for example). I have also suggested as the loop BW gets to large we start introducing other noise sources and lose signal by tracking out the carrier such that an optimization in loop BW is involved. The OP is making the good point that it would be impossible to track out the carrier if we are only sampling 1 sample/symbol. So the question is what others factors would limit the maximum carrier tracking loop bandwidth, or would the go be to provide the maximum bandwidth that can be achieved in a loop with an update rate $R$ of a sample/symbol (i.e. < $R/2$).

Below I show more detailed simulations of this type of decision directed carrier tracking loop. Specifically we see the closed loop transfer of noise from the input to just prior to the decision as a high pass filter (as the loop tracks out the slower changes in phase as part of tracking the carrier). We note that with a narrow loop bandwidth the gain levels out close to 0 dB at the high end of the Nyquist bandwidth, while in the wide loop bandwidth the gain of the higher frequencies is +6.8 dB, which leads to noise enhancement.

The expected input noise will typically rise significantly as dominant phase noise as we approach the carrier (low frequency offsets) as is typical for local oscillator phase noise. However at some offset from the carrier the decreasing phase noise will intersect the noise floor which will then dominate from other white noise sources such as amplified thermal noise and quantization noise in fixed point implementations. Thus we see how a trade can exist with the loop bandwidth setting where the savings from further attenuation of lower frequency noise is offset by the noise enhancement from the wider loop bandwidth setting. Note that the decision directed phase detector is equally sensitive to AM and PM components on the signal, and the AM noise components of the input signal will be translated to PM components on the NCO in the loop and will thus not cancel out the AM noise but add to it as an uncorrelated noise source. Ultimately if an optimization is desired, and to confirm if this particular noise contribution will even affect the receiver,the simulation and analysis would need to be done with the expected noise profile for the specific system being analyzed, and importantly while operating at senstivity. So there is not a one answer fits all but hopefully this makes it clear how the SNR could be degraded by the carrier tracking loop if the bandwidth is set too high and what the considerations are for analyzing this.

Simulation Result with Phase and AWG Noise

Here is a complete simulation of the Decision Directed Carrier Tracking Loop with a 2nd Order PI filter (as a typical implementation) showing how the loop bandwidth can be optimized for pre-decision SNR performance, specifically showing conditions where the noise can degrade as the loop bandwidth continues to get wider, even if it is not at the point of tracking out the modulation.

Given the following Carrier Tracking Loop (Applicable to BPSK, QPSK, M-PSK and QAM modulations):

The following is the Loop Model and equation for the open loop gain:

$P$: Proportional gain constant

$I$: Integral gain constant (per sample with time in seconds)

$IT$: Integral gain consant (normalized with time units in samples)

$A$: complex magnitude of samples at discriminator input

$K_d$: Discriminator gain includes parasitic 1 sample delay in loop

$H_{LF}(z)$: PI Loop Filter

$K_{NCO}$: NCO gain, control word given in frequency of Hz (this way we can monitor frequency directly)

$$H_{LF}(z) = \frac{Pz + IT-P}{z-1}$$

$$K_{NCO} = 2\pi\frac{z}{z-1}$$

$$K_d = \frac{A^2}{z}$$

$$G_{OL}(z) = 2\pi PA^2\frac{z-\bigg(\frac{P - IT}{P}  \bigg)}{(z-1)^2}$$

With this Loop Model I created three test cases, for all cases $A=1$ and $T=1$ (normalized magnitude and time is in samples):
$$\begin{bmatrix}  \text{Loop BW} & -3 \text{ dB BW}       & \text{Settles In}  & I     & P    & G_{CL}(z=\pi)   \\
\text{---------}&\text{---------------}&\text{-------------}&\text{-------}&\text{------}&\text{----------------}\\
\text{Narrow}  & 0.011 \text{ cyc/samp} & 200 \text{ samples} & 1e-4  & 7e-3 & 1.045 \\
\text{Medium}  & 0.056 \text{ cyc/samp} & 20 \text{ samples} & 0.01  & 0.07 & 1.257 \\
\text{Wide}  & 0.202 \text{ cyc/samp} & 5 \text{ samples} & 0.072  & 0.21 & 2.206 \\
 \end{bmatrix}$$

I then created a noise sample of $2^{15}$ complex samples with the following phase noise target values for the PSD:

Freq $1=2\pi$,   Phase Noise (dBc):

1e-2 cycles/sample, -15 dBc

1e-1 cycles/sample, -45 dBc

0.5 cylces/sample, -59 dBc

To this phase noise data I added AGWN with a total power of -40 dBc. This created experimental data of complex samples with both AM and PM noise components, with the phase noise contributing AM only with increasing density toward the lower frequencies, and AWGN contributing AM and PM components equally.

The loop performance was characterized by determining the closed loop transfer function from the input of the phase rotator to the output of the phase rotator which is given by:

$$G_{CL}(z) = \frac{1}{1+G_{OL}(z)} $$

It is clear that signal component, here normalized to 1 would not be affected in magnitude by the phase rotator, so we can assess SNR from the noise directly after it has passed through the above transfer function.

The Decision Directed Phase Detector responds equally to small scale AM and PM changes (A small change in amplitude can not be distinguished from a small change in phase), so any AM components would get translated to (uncorrelated) PM noise at the output of the phase rotator as the loop tries to correct phase offsets that don't exist. Therefore the resulting total noise at the output taken at the phase rotator output would be:

$$S_o(\omega) = S_{am}(\omega) + |G_{CL}(\omega)|^2S_{\phi}(\omega) + |G_{CL}(\omega)|^2S_{am}(\omega) $$

Where $S_o(\omega)$ is the power spectral density of the noise alone. Which means all the AM noise original on the signal would still exist without modification since the loop does not correct for any AM, the PM noise would be filtered by the loop, but the AM noise would also be translated to PM by the loop and this component of translated AM noise would also be filtered by the loop.

So accounting for all that, using this experimental data, here were the results after passing the data set through the narrow, medium and wide bandwidth loop filter responses:

$$\begin{bmatrix}  \text{Loop BW} & \text{Total Integrated Noise}  \\
\text{---------}&\text{-----------------------------}\\
\text{Narrow}  & -64.6 \text{ dBc} \\
\text{Medium}  & -69.8 \text{ dBc} \\
\text{Wide}    & -67.1 \text{ dBc}\\
 \end{bmatrix}$$

I assumed the optimum loop bandwidth for minimizing noise would be around the corner frequency where the phase noise intersects with the white noise floor but haven't confirmed this, nor did I do a detail sweep beyond these three test cases. I wanted to show an example that demonstrates how wider loop bandwidths can cause the noise level to start increasing as bandwidth gets wider and for this loop it is because of noise enhancement for wider loop bandwidths. We can see how this must be the case given the two closed loop zeros at $z=1$, and that the poles move closer to $z=-1$ as loop bandwidth increases.

Also to note, the resulting overall noise levels in this test case are inconsequential for demodulating most communication waveforms, so it would make no difference to degrade the tracking noise by 3 dB so that we could operate at the wider loop bandwidth, with all its other benefits. The point here was to show that it can be optimized but also that it very much depends on the overall system and there is not a one answer for all other than to evaluate this on a candidate system with similar detail.

The above is the concluding data and remarks, but further plots are shown below with regards to the test data, loop step responses and closed loop transfer functions.

Experimental Data - Phase Noise Only (time domain)

Power Spectral Density of Phase Noise Data

Power Spectral Density of Phase Noise with AWGN Data

Closed Loop Step and Frequency Response From Input to Phase Rotator Output


Filter Cascade: Additions and Multiplications per input sample Recently I was reading "Optimum FIR Digital Filter Implementations for Decimation, Interpolation, and Narrow-Band Filtering) by Crochiere and Rabiner (1975).
To estimate the number of decimation Stages and factors, they are targeting to optimize the total number of multiplies and adds (per second) calling this number $R_T$. With $R_T$ being sum of multiplies and adds of each stage:
$R_T = \sum_i^K R_i$.
When playing around with Matlab, I found the Matlab function designMultistageDecimator, in one  example they compare a multistage approach with a single filter by calculating the cost (with the function cost ), here however, it is multiplication and addition per input sample.
I would have assumed that the cost behaves in the manner, that it is the total sum of all cost per stage.
For the single filter, which will be then compared to a cascade, I understand the result, having 753 coefficients and a decimation factor of 48, that roughly there are 16 mults and 16 adds per input sample.
But when it comes to the filter cascade, in the example given by Matlab

Stage one has 17 coefficients (13 non zero) and a decimation factor of 3
Stage two has 11 coefficients (7 non zero) and a decimation factor of 2
Stage three has 15 coefficients (9 non zero) and a decimation factor of 2
Stage four has 63 coefficients (49 non zero) and a decimation factor of 4

The estimated multiplications are ~7.2 and additions ~6.6 per input sample.
This confuses me, I would have assumed that is roughly ~25 multiplication and ~25 additions. Just taking the last stage having more or less 50 non zero elements and a decimation factor of 4 would add at least 12.5 operations, or not?
So my questions are,

How are the costs, in the sense of multiplies and adds computed for a cascaded design in matlab using the function cost? (Thanks @TimWescott)

Why is it not the sum of multiplication and add per stage? A sample at the output has to see at least on o the polyphaser trees of the last stage having at least 12 coefficients or not?

Follow up question:

After the answer by @TimWescott, I understand now the math behind the result of the cost funciton. However, I am asking myself how to interpret the output of the cost function? Is on sample, given this example, really only seeing ~7 multiplications / ~7 additions?
In therms of dsp-terminolgy what would the calculation of operations int the following way:

13/3 + 7/2 +9/2 +49/4 = ~25
represent?

Their numbers work out right if you do the stated decimation before the stage, and you pay attention to the fact that decimations are cumulative.
So stage one operates at 1/3 of the input rate, stage two operates at 1/6, stage three operates at 1/12, and stage four operates at 1/48.  If you scale the adds and multiplies by those numbers, then their stated numbers work.


What is the relation between the PSDs of filter input and output called? $R_Y = |H|^2R_X$ If a wide-sense stationary signal $X$ is fed to an LTI filter with the transfer function $H$, the power spectral density (PSD) of the output $Y$ can be expressed as:

$$R_Y(f) = \left|H(f)\right|^2R_X(f)$$

where $R_X$ denotes the PSD of $X$.

Does this relation have a common name?

I don't know the name of the relationship, but $\vert H(f)\vert^2$ is
called the power transfer function of the LTI system.  The output
power spectrum is the input power spectrum multiplied by the
power transfer function, just as for deterministic signals, the output
spectrum is the input spectrum multiplied by the transfer function
$H(f)$.


Removing high pitch noise embedded in audio spectrum I have an audio clip. Here is the clip in wave format

https://www.dropbox.com/s/542jgpfk7bfkt01/sample.wav?dl=0

I chose the following area of time domain signal

The spectrum of this signal shows maximum content at 2371 Hz which is inside voice band.

I applied notch on this frequency but i can still hear this noise.

Can any body tell how to remove it completely so that it is not hearable?
same is what i want with noise in beginning of clip.

It will be bonus for me if someone can tell why this high pitch noise occur in audio?
I will be interested in knowing DSP techniques which I can apply in MATLAB rather than sound processing software.

Can any body tell how to remove it completely so that it is not hearable?

As suggested already, the signal has harmonics.

You can diminish it further by keep filtering the rest of the peaks of the spectrum with a notch filter, or using the noise reduction option in Audacity.

The latter option is better when the interference is over a useful part of the recording. To apply that filter, select a part of the recording that is representative of the noise profile you are dealing with. Audacity will use this to extract the noise spectrum and then create a filter that cancels it. You can then take that profile and apply it to other points in the recording to reduce its impact.

  It will be bonus for me if someone can tell why this high pitch noise occur in audio?

It sounds like interference and it is likely that it is picked up as a square wave. It definitely doesn't sound like a hum, more like a high pitch pipe organ sound with two or three components.

It might be coming from unshielded cables, strong sources nearby an amplifier (e.g. WiFi router / dongle, LED lightstrip, TV, other) or if it is an audio USB dongle, of dubious built quality, it might be picking up interference from the data exchange going on on the bus. It's hard to tell without knowing what is this a recording of (?).

Hope this helps.


Power spectral density of sinusoidal signal in noise In a literature I face with this input and power spectral density (PSD)
$$x(t)=s(t)+n(t)=A\cos\left(\omega_c t +\phi\right) + n(t)$$

first I want to know

How can we find PSD of $\cos\left(\omega_c t +\phi\right)$?(I have found the autocorrelation function of $\cos\left(\omega_c t +\phi\right)$ but when I apply the PSD with Fourier transform it doesn't return the same result as the figure)
Why do we divided PSD of $\cos\left(\omega_c t +\phi\right)$ by $R_L$?
What is a bandpass filter's transfer function in time and frequency? And how can we find PSD of a bandpass?

Sorry for the simple question but I think I got confusing.

It looks like you are looking at the result of a combination of things.

First of all let us observe the image and try to understand it. The shredded horizontal line with value $\frac{N_{in}}{2B}$ represents $n(t)$ in the frequency domain and the signal $s(t)$ is represented by the two vertical arrows valued $\frac{A^2}{4R_L}$. Two arrows one in the sinus's frequency value and its mirror in the negative side is a representation of a sinuous wave in the frequency domain. The cosinuous wave will have the a similar representation as the sinuous wave in the frequency domain.I will reffer you to the table here

A BPF is a band pass filter and it is a combination of a low pass filter (will only pass frequencies up to a specified $f_0$ frequency) and a high pass filter (will only pass frequencies up from a specified $f_0$ frequency). In short, a BPF will only pass frequencies from a specific band bounded between some $f_1$ and $f_2$ as in $f_1\leq f\leq f_2$. The BPF is represented as a window in the frequency domain. You can see in your image the window with $B$ width around the sine arrows. This is also the reason the noise $n(t)$ is shredded, The filter is cutting it off outside the window's boundary. In the time domain a BPF will be represented by a $sinc()$ function, as can be seen in the table in this page.

I think your question needs additional data, since it lacks an explanation of what exactly $R_L$ is. I am assuming it is a resistor and this is some sort of an electronic circuit but if you want any further help you need to supply more detailed on the system.


Fixing edge/pixel defects in image I have a large number of images that are supposed to look something like the top part of the image:

Instead, some have defects in them, as in the bottom part of the image. As you can see, there are 2 types of defects:

Outliers: These are pixels (or small areas) of one colour located in a vast area of another colour.
Missing edges: These are black areas on the edges between two (non-black) colours.

There's also a third defect in some images, which is blur on the edges, but I think I can deal with that myself.

I'm looking for a way to automatically fix these images so that:

The outliers get replaced with the colour they're located in (so a red pixel in the green area gets painted green).
The missing edges get coloured in with the outer colour (so the black edge on the boundary between green and red gets painted green).
The rest of the image is left as is.
No additional defects (such as blur or image noise) are introduced in the process. The only colours present in the 3 ellipses should be (0,255,0), (255,0,0), (0,0,255).

Is such an automatic process possible using some kind of filter-based method or would I have to do this with some clever if-then rules for handling pixels based on their neighbourhoods?

What you are looking at for this problem is a combination of "Image Segmentation" and "Spatial Analysis".

The image segmentation part of this particular task is very simple because each pixel "class" is already labeled by a colour. You might need an additional labeling step because (for example) the blue class contains two disjoint areas, but this is easily dealt with via a simple connected component labeling. (We will come back to this in more detail).

The more challenging bit here is what corresponds to relationships expressed in natural language as "the blue areas inside the red areas", or "the black areas adjacent to red areas" and so on. These are relationships that correspond to mathematical "queries" such as "Point on line", "Point within area", "Area within area" and so on.

For example, to determine if a point is within a circle, all you have to test for is whether the distance from the centre of that circle to the point in question is less than the radius of the circle. To test a point within a rectangle you would test if its coordinates are within the boundaries of the rectangle and so on.

As you may have guessed, this does not sound easy to generalise (e.g. any kind of polygon, whether in 2D or 3D or higher dimensions etc). Fortuneatly, there are a lot of good libraries out there that implement geometrical primitives, operations and relationships for us to make such tasks easier.

In the case of Python there is Shapely. Shapely contains primitives such as Point, Line, Polygon and others as well as operators and relationships you can use to "query" their spatial relationships.

Here is a very simple example that answers part of your task: Is one area within another area?
Notice here, I simply wrote "area" as this is defined by an "unordered" set of pixels similar to what you would get out of simple queries such as "Give me all pixels in an image that are red".

    import shapely.geometry

    Q = shapely.geometry.MultiPoint([(u,v) for u in range(0,10) for v in range(0,10)])
    W = shapely.geometry.MultiPoint([(u,v) for u in range(5,8) for v in range(5,8)])

With these two lines we have created two MultiPoint (areas composed of many Points) "areas": One big one from (0,0) to (9,9) and one small one from (5,5) to (7,7). Furthermore, the small "patch" is contained by the bigger "patch".

Let's probe that relationship

    Q.contains(W)

This should return True

   W.contains(Q)

This should return False.

You might be wondering what would happen if what we called the "small patch" contained two disjoint patches. In that case, Shapely would respond with False if only one of the patches was found within a testing area.

The black areas are a bit more challenging. You might think that their case is covered by the .touches() relationship but "touching" implies some sort of overlap and there are no criteria here by which to make those areas somehow to overlap. But what you can probe instead is distance.

In the case of MultiPoint objects, the distance between two objects is the shortest possible distance between their points. Again, here is a minimal example:

    import shapely.geometry

    Q = shapely.geometry.MultiPoint([(u,v) for u in range(0,10) for v in range(0,10)])
    W = shapely.geometry.MultiPoint([(u,v) for u in range(0,10) for v in range(12,22)])

Notice the differences in the way the "boxes" are setup here. To get the distance between them:

    Q.distance(W)

This should return 3.

You can see here that the black areas distance to blue, red and green is 1 (they are adjacent). Therefore, to locate which parts is a given black part adjacent to, you just assess its distance from the parts of interest.

So, with this bit out of the way, let's get back to how you would segment/label your pixel masks.

If your regions where not overlapping it would be extremely simple to extract the pixel coordinates that correspond to all pixels that have a specific color and create MultiPoint instances from them. Here is an example using Python's numpy and Shapely.

    import skimage.io
    import numpy
    import shapely.geometry

    # Read the image
    I = skimage.io.imread("someRRGBimage.png")
    # Exrtact the pixels that are equal to a specific value
    pixel_bunch = numpy.where((I[:,:,0]==the_red) & (I[:,:,1]==the_green) & (I[:,:,2]==the_blue)) # Notice here the_red, the_green, the_blue are the components of the area's colour.
    # Create a MultiPoint patch
    point_patch = shapely.geometry.MultiPoint(list(zip(pixel_bunch[0], pixel_bunch[1])));

But if you try to select the blues or the reds out of the image you are providing, you will get one large area, plus disjoint little ones and if you try to pass those through the spatial relationship analysis we are talking about above, they would fail.

The simplest way to split them apart is to use connectivity labeling. This works by scanning a binary image and assigning a different number to regions that are not connected. This basically inserts an extra step to the extraction of the pixel values above, as you would first have to create the mask image for a specific colour, then label that mask and finally create the regions based on the identified (and by now isolated) regions.

It's far easier to demonstrate this over grayscale images:

    import skimage.io
    import skimage.measure
    import numpy
    import Shapely.geometry

    # Read the image
    I = skimage.io.imread("someRRGBimage.png")
    # Work on the Red channel
    mask = (I[:,:,0]==255)
    labeled_mask = skimage.measure.label(mask)
    pixel_bunch = numpy.where(labeled_mask==some_patch_id) # Notice here some_patch_id is the zero based index of some isolated patch that has been labeled in the image.
    # Create a MultiPoint patch
    ...
    ...

It is not a difficult task but it sure is messy.

Finally, two notes:

You will probably have to re-colour those areas between the red and the blue that are currently black and leave black for the "background".
The complexity of deriving those relationships scales with the size of the regions, sometimes taking tens of seconds (or more) for very large MultiPoint sets. In that case, it might be more practical to work with polygons. This inserts yet another step just before deriving point_patch where you first derive the boundary of the identified pixel area and then instantiate a Polygon from that boundary. Operations between polygons will be much quicker.

Hope this helps.


Approximating pixel location in scale space I'm implementing a scale-space with o octaves and s scale levels in each octave.
Each octave is half the size the previous.
I have keypoints found in each octave, I need to approximate the real location of the keypoint to their location in the original image (before the reduction). I need it to be more accurate then just doubling the coordinates.
I know how to approximate sub-pixel in each octave, but I need "between" the octaves.

Thanks
matlabit

First, you need to ensure that the pixel has minimal value in its 8-neighborhood. SIFT, SURF and other keypoint detectors filter these pixels in step called "non-maximal suppression". This is basically a necessary condition for second-order approximation we will use to determine sub-pixel location.

If you imagine the scale space slice $I$ as a 3D surface, the pixel lays in a valley sorrounded by pixels with higer value. The shape can be well approximated by a parabola surface and the corrected pixel position is at its minimum.

The pixel value $I(p)$ at location $p=(x,y)^{\textbf{T}}$ can be approximated by second order Taylor expansion:

$$I(p+\delta)\approx I(p) + \delta^{\textbf{T}}g+\frac{1}{2}\delta^{\textbf{T}}H\delta$$

where $g$ and $H$ are gradient (2-vector of first order derivatives) and Hessian ($2\times 2$ matrix of second-order derivatives) at $(x,y)$.

The $\delta$ is a sub-pixel correction factor we want to compute.

Since the parabola has minimum at point where its derivatives are equal to zero, we obtain derivative of right-hand side of the above Taylor expansion (with respect to $\delta$) and place it equal to zero:

$$g+H\delta=0$$

From this we can compute the correction vector:

$$\delta=-H^{-1}g$$

Finally, the desired sub-pixel location is given by $p + \delta$.

The gradient and Hessian are obtained using finite differencing:

$$g_{1}=\left( I(x+1,y)-I(x-1,y) \right)\cdot 1/2$$
$$g_{2}=\left( I(x,y+1)-I(x,y-1) \right)\cdot 1/2$$
$$H_{1,1}=I(x+1,y)+I(x-1,y)-2I(x,y)$$
$$H_{1,2}=H_{2,1}=\\ \left(I(x+1,y+1)+I(x-1,y-1)+I(x+1,y-1)+I(x-1,y+1)\right)\cdot 1/4$$

So the algorithm is:

Check if the pixel is a local minimum (sorrounded by pixels with higher value).
Compute $g$ and $H$ using finite differencing formulas above.
Compute $\delta$.

You can repeat steps 2. and 3. for extra accuracy, but just one step is fine for our needs.

The value at sub-pixel location $I(p+\delta)$ can be obtained by well known bilinear or bicubic interpolation.

We did the interpolation with respect to $x,y$, but you can perform this interpolation with respect to scale as well. This will lead to 3-vector gradient and $3\times 3$ Hessian, but the Taylor expansion formula holds.

You can take a look on OpenSURF implementation in C++ or C# which has this implemented.

The above derivation is actually a Newton's Method used in multivariate optimization.


Peaking filter (loudness) cutoff frequency in an equalizer An equalizer we have consists of Loudness ,Bass ,Mid and Treble .The loudness is nothing but a Peaking filter with cutoff frequency 20Hz .I read this in a book that outlines the EQ design.
I was wondering why the Loudness block which is basically a Peaking Filter should have a cut off frequency of 20Hz ?? Any intuition behind this

The purpose of the loudness filter is to compensate for the level dependency of the a human's "frequency response". Most music is mixed at pretty high levels, say at 80-90 dB SPL. If you play this back at much lower levels, say 40 dB SPL, it will sound very bass deficient. The loudness filter compensates for this by boosting the bass at lower levels.
In a "good" implementation the loudness adjusts automatically based playback volume setting and content energy: the lower the playback volume the more bass boost is needed to maintain the same perceived spectral balance between low and mid frequencies.
The shape of the loudness filter can be derived from the Equal Loudness contour of humans (see https://en.wikipedia.org/wiki/Equal-loudness_contour). In essence the loudness filter is the difference between the equal loudness curve at the playback level and the equal loudness curve at some reference level, all normalized to 0 dB at 1 kHz (or something like this).
The target filters have a continuously rising bass response. That's not easy to implement with "standard" EQ blocks. A very low frequency peak with the right choice of center frequency and Q can be a fairly good match in the region you care about. You use a 20 Hz peak filter because it has the right shape from 40Hz-100Hz and you don't care much what happens at 20 Hz since there almost no content there.


Nyquist Plot for transfer functions with poles at the origin I'm learning Nyquist plots and something has been seriously bugging me when treating poles or zeros in the origin. Nyquist plots obtains information based on the argument principle which states

"If f(z) is a meromorphic function inside and on some closed contour C, and f has no zeros or poles on C, then

$${\displaystyle \oint _{C}{f'(z) \over f(z)}\,dz=2\pi i(N-P)} \oint_{C} {f'(z) \over f(z)}\, dz=2\pi i (N-P)$$
where N and P denote respectively the number of zeros and poles of f(z) inside the contour C, with each zero and pole counted as many times as its multiplicity and order, respectively, indicate."

So we just ignore the fact the transfer function contains a pole over the contour $C$? How can we treat them normally when it clearly violates the argument principle?

We don't ignore poles on the contour. As mentioned in a comment, poles are avoided by modifying the contour as shown in the figure below, where a contour appropriate for a pole at $s=0$ is shown.

Fig. 1: Nyquist contour for a pole at $s=0$ (from "Modern Control Engineering" by K. Ogata).

The contour moves around the pole along a semi-circle centered at the location of the pole. The radius of that semi-circle approaches zero, such that the whole right half-plane is enclosed by the resulting contour. Note that by choosing the contour in this way, a pole on the imaginary axis is outside the contour, and will not add to the encirclements of the origin in the Nyquist plot.

Of course we could also move along a semi-circle in the left half-plane to avoid a pole on the imaginary axis. In that case, the pole would be inside the contour.

As an example, consider the function

$$F(s)=\frac{(s+2)^2}{s(s+1)}$$

It has a double zero in the left half-plane, one pole in the left half-plane, and one pole on the imaginary axis at $s=0$. If we use the contour shown in Fig.1, we get the Nyquist plot shown in the right-hand side figure below (the corresponding contour is shown on the left).

There is no encirclement of the origin, in agreement with the fact that there are no poles and zeros inside the contour. Note that due to our choice of the contour, the pole at $s=0$ is outside the contour.

If we choose a different contour with a small semi-circle in the left half-plane to avoid the pole at $s=0$ (left-hand side figure below), the pole at $s=0$ is inside the contour, and, consequently, the Nyquist plot shows one counter-clockwise encirclement of the origin (right-hand side figure below), corresponding to one pole and no zeros inside the contour.

In sum, poles on the imaginary axis are avoided by moving along semi-circles of infinitesimal radius, and, depending on whether the semi-circle is in the right or the left half-plane, the poles on the imaginary axis are either outside or inside the contour, which is reflected in the Nyquist plot by the number of encirclements of the origin.


Prepare Buffer for Inverse FFT How should buffer value at N/2 get initialized prior to performing an inverse FFT on a buffer, size N, that contains a complex frequency response in only values 0 to N/2-1?  It is pretty clear the buffer needs to be made symmetrical for prior to calculating the inverse FFT per table 12-1 of http://www.dspguide.com/ch12/1.htm.  What is not clear is how to initialize the value at N/2.  Is it related to value at 0?

Normally a buffer of N points of data in the frequency domain covers, much like a two's complement number, -N/2 to +N/2-1.  So when you say that it doesn't cover N/2, it probably doesn't matter.  If it does matter and there is a "place" for N/2, just do as Paul suggested and set it to 0.


Image Pyramid Without Decimation I have an image processing algorithm which uses Gaussian and Laplacian Pyramid Decomposition.
At each level of the pyramid the algorithm decimates the previous level by a factor of 2.
The LPF the algorithm uses is Gaussian.

I was wondering whether I could get the same results without the decimation step (By altering the LPF).

As far as I know, the decimation step is equivalent of cutting by half the bandwidth of the LPF.

What I tried, with no success, is to double the STD of the LPF at each step and skip the decimation.

Is there another / better way to do it?
Thanks.

You can generate a scale volume from an image by convolving it with Gaussians of increasing $\sigma$. Or, equivalently, by convlolving it repeatedly with the same Gaussian.  It just won't be a pyramid. :)

One reason to generate a pyramid by downsampling is to reduce the amount of computation, because once you apply a low-pass filter the full resolution becomes redundant.

Also, are you computing the Laplacian pyramid by using the difference of Gaussians?  In that case, you need to downsample each subsequent scale level, and then upsample it before you take the difference with the previous level.


Determine g(t) for minimal symbol error probability, given h(t), K=1 and Minimum distance receiver Determine $g(t)$ for minimal symbol error probability, given $h(t)$, $K=1$ and Minimum distance receiver.
$h(t)=   \begin{cases}
      0, & \text{$t<0$}\\
      2, & \text{$t\in(0\;, 10^{-3})$}\\
      1, & \text{$t\in(10^{-3}\;,2\cdot 10^{-3})$}\\
      0, & \text{$t>2\cdot 10^{-3}$}
    \end{cases}  $
This is ISI-free if $h(t)*g(t)|_{t=kT}=
    \begin{cases}
      1, & \text{if $k=0$}\\
      0, & \text{if $k\neq0$}
    \end{cases}  $
Using a matched filter as $g(t)=ch^*(-t)$ the convolution becomes
$$(g*h)(t)=\int_{-\infty}^\infty{h(\tau)g(-\tau)d\tau }=c\int_{-\infty}^\infty h(\tau)h^*(\tau)d\tau=1$$
How do I calculate the $\int_{-\infty}^\infty h(\tau)h^*(\tau)d\tau=1$ integral or the pulse energy?
From what I understand it is the same as the pulse energy (which I don't know how to calculate either).
From here I understand how to determine $g(t)$ and it will be the $c$-scaled and time inversed variant of $h(t)$, but I don't really understand how to calculate the pulse energy.
Is $E_h=\int_{-\infty}^{\infty}|h(\tau)|^2d\tau$$\rightarrow$
$E_h=\int_{0}^{10^{-3}}2^2d\tau+\int_{10^{-3}}^{2*10^{-3}}1^2d\tau$ the correct way?

Your computation using the energy of $h(t)$, $E_h$, to find the value of $c$, is correct. But note that you are only forcing the zero ISI condition for $k=0$. You still need to check the zero ISI criterion whenever $k\neq0$, but there is no mention of the value of $T$, and we can't derive what happens in those cases until we know the value of $T$.
For a general approach, I'd recommend going to the basics and not taking shortcuts. One key problem in your definition is that your convolution integral is not taking into account the time instant $t$, and therefore, you only compute one value of the convolution (for $t=0$), instead of all the values that are relevant (namely, $kT$ for all $k$).

The convolution of $h(t)$ and $g(t)$ is defined as
$$p(t) = (g* h)(t) = \int_{-\infty}^{\infty}g(\tau)h(t-\tau)\,d\tau = \int_{-\infty}^{\infty}h(\tau)g(t-\tau)\,d\tau$$
and we need to satisfy the Nyquist first criterion for zero ISI,
$$p(t)\Big|_{t=kT}=(g* h)(kT)=\left\{\begin{array}{rl}1,&k=0\\0,&k\neq0\end{array}\right.$$
Lastly, the matched filter is defined as $g(t)=ch^*(T-t)$ in its causal version, and $g(t)=ch^*(-t)$ in its non-causal version. For simplicity, let's use the noncausal version:
$$\begin{align}
p(t)& = (g* h)(t) = c\int_{-\infty}^{\infty}h(\tau)h^*(\tau-t)\,d\tau\\
\end{align}$$
In the particular case of $t=0$ (equivalently, $k=0$), we have
$$\begin{align}
p(0)& = (g* h)(0) = 1\\
&=c\int_{-\infty}^{\infty}h(\tau)h^*(\tau)\,d\tau= c\int_{-\infty}^{\infty}|h(\tau)|^2\,d\tau\\
&= c\Bigg(\int_{0}^{10^{-3}}|2|^2\,d\tau+\int_{10^{-3}}^{2\cdot10^{-3}}|1|^2\,d\tau\Bigg)=c\big(4\cdot10^{-3}+10^{-3}\big)=5c\cdot10^{-3}\\
&\Rightarrow c=10^3/5 = 200
\end{align}$$
Now, for verifying that the system satisfies the zero ISI criterion, we need the value of $T$. But you can follow a similar procedure.


FFT resolution and sampling frequency I am computing FFT of the data which I am receiving from an ADC with SPI protocol. Since resolution of FFT is related to $F_s/N$ where $N$ is the FFT length and $F_s$ is the sampling frequency of ADC, when I lower sampling frequency, it should decrease frequency range that I see due to Nyquist–Shannon sampling theorem and conversely it should also increase my FFT resolution.

I am plotting this FFT in realtime and I've given 1 Hz sinusoidal signal via arbitrary signal generator and tested various sampling frequencies to see which point 1 Hz equals to.

╔═══╦═══════════════╦═══════════════╗
║   ║ Sampling Freq ║ Point of 1 Hz ║  0. point included
╠═══╬═══════════════╬═══════════════╣
║ 1 ║ 20k           ║ 3             ║
║ 2 ║ 10k           ║ 6             ║
║ 3 ║ 1k            ║ 9             ║
║ 4 ║ 100           ║ 9             ║
╚═══╩═══════════════╩═══════════════╝

No matter how I lower sampling frequency, it doesn't go further than 9. By the way, my FFT length is 1024.

One issue is that, for the first two sampling frequencies, you are not sampling a full wave of the signal.

Below is a plot of the time domain signal (1024 samples) and the resulting FFT.  For the first two, the DC bin (0 index) will capture all the energy.  For the 1kHz example, the first bin (the one after 0 index) will capture it.  Bin 10 is where the peak is for 100 Hz sampling rate.

Perhaps your sampling isn't happening quite the way you expect?

R Code Below

# 43563

N <- 1024

fs <- c(20000, 10000, 1000, 100)
res <- c(0,0,0,0)
idx <- 1

par(mfrow=c(4,2))
for (freq in fs)
{
  res[idx] <- freq / N
  idx <- idx + 1
  t <- seq(0,N-1) / freq
  x <- sin(2*pi*t + 0.892340923)
  plot(x, type='l')
  title(paste('time domain for ', freq))
  x_fft <- abs(fft(x))
  plot(x_fft[1:100], type='l')
  title(paste('100 bins in frequency domain for ', freq))
}


What is the most efficient type of software bandpass filter? I'm writing a C++ application to extract four modulated carrier signals at four separate frequencies, from a single input signal. I've designed four FIR filters at t-filter.appspot.com but they run very slowly, even when I reduce # of taps to 80.

Is there a better option? People seem to point to biquad filters – are these different from FIR filters, and do they run more quickly in software?

Frequency domain filtering (FFT), as suggested by some comments, is definitely wrong -- it's even slower, or same speed at best! A recursive filter (IIR) is the fastes possible solution. If you choose a typical second order filter (called biquad in engineering slang) of Butterworth type and do your math right (factoring out coefficients) you only have 3 multiplications and five addition.

Edit: link to classics page: filter cookbook


what is endpoint in image processing? what is endpoint? I have a sample code to detect endpoints . it finds endpoint with sample, but when i run for my sample tif image, it does not well. so I want know the exactly meaning of endpoint in image processing.

This concept is related to the topological skeletonization of an object in shape analysis. It produces a thin sert of lines equidistant to the boundaries, preserving  geometrical and topological properties of the object. Some use the term medial axis. Some points possess a specific interest: branch or junction points (where sereval lines meet, i.e. a point with more than two neighbors) and termination, terminating or end points (a point with less than two neighbors), located at extremities.


Sampling orthogonal signals Suppose I have two signals $f,g$ which are orthogonal in the sense of the usual inner product on $\mathbb R^d$. Then if I'm not mistaken these signals have disjoint spectral support, i.e their fourier transforms have disjoint support. Suppose further that I know the spectral supports aswell. How can I make use of this information in order to sample only one of the signals, i.e, sample so I don't "see" the other at all?

The trivial solution seems involve filtering out the irrelevant frequencies using an ideal filter, but that's not an option; I am restricted to "physical" solutions.

It's fine to assume $f,g$ are well-behaved.

I'm afraid that your premise is wrong. The Fourier transforms of two orthogonal signals do not in general have disjoint support. You know that two signals $u(t)$ and $v(t)$ are orthogonal if

$$\int_{-\infty}^{\infty}u(t)v^*(t)dt=\int_{-\infty}^{\infty}U(f)V^*(f)df=0\tag{1}$$

where $U(f)$ and $V(f)$ are the Fourier transforms of $u(t)$ and $v(t)$, respectively, and $*$ denotes complex conjugation. One possibility for (1) to be true is indeed that the product $U(f)V^*(f)=0$, but this is of course a sufficient condition, not a necessary one.


Finding the fundamental period by using four elementary operations (+,-,x,/) The idea is to find the fundamental period of any periodic and bounded signal. This is all that we know about the signal. Nothing else. Only allowed operations are: $+,-,\times , \div$. We can also access any value of the signal that we desire.

How can a method be developed which finds or estimates the fundamental period of the given signal?

Thank you in advance.

Well eventually operations such as the Fourier Transform and autocorrelation are defined with those four operations.

Are comparisons and conditional decision-making allowed?  Are you allowed to pick a pair of numbers from an array based on one of those numbers being a maximum?

My other answer is a place to begin.


Principal Component Analysis as a step between STFT and MFCC For computing Mel Frequency Cepstral Coefficients you can use already calculated STFTs as a basis and perform the Mel frequency mapping on it.

My question: Does it make a significant difference if I calculate the STFTs, perform a PCA-transform on them and then calculate the MFCC compared to computing the PCA-transform at the very end on the MFCCs?

In both cases we have two dimensionality reductions (PCA and MFCC) and I'm not sure weather the order makes a difference.

One dimensionality reduction method (PCA) is linear, the other (MFCC) is  not. It  is thus unlikely that there is no difference at all since they probably do not commute. One serves more decorrelation, the other feature extraction.

The importance of the order thus depends a lot on your signal class, and what your purpose is.


Daubechies Wavelet and Matlab I want to use waverec to evaluate a linear combination of the scaled and shifted wavelet  of the form

$$ \sum_{i=1}^n\sum_{j=0}^{2^n-1}d_{nj} 2^{-n}\psi_{ij} (\frac{t-j}{2^n}) +c_{00}\phi(t)  $$
at $t=k2^{-n} $ for $k=0\dots2^{n}-1$. I know how do it for the Haar wavelet, namely

[c,l]=wavedec(zeros(1,8),3,'haar')
waverec(2^(3/2)*[c00 d10 d20 d21 d30 d31 d32 d33],l,'haar')

But how can you do the same thing for Daubechies wavelets? The problem being that the array c does not have 8 entries but 28:

[c,l]=wavedec(zeros(1,8),3,'db4')

c =

  Columns 1 through 14

     0     0     0     0     0     0     0     0     0     0     0     0     0     0

  Columns 15 through 28

     0     0     0     0     0     0     0     0     0     0     0     0     0     0

l =

     7     7     7     7     8

The Question is related to Number of Daubechies coefficients, but I do not understand the answer.

The wavdec function does not produce a DWT which is non-redundant. Instead, wavdec takes you N levels down the decomposition process. In other words, you shouldn't expect to get 8 coefficients back.


Fourier transformed acceleration from Fourier Transform of velocity and position using Omega arithmetic Let's say the position of an object is given by simple sine function. By elementary calculus, I can calculate the acceleration in the time domain and find its Fourier transform.
I can also calculate the Fourier transform of velocity and multiply it by frequency or multiply the Fourier transform of position by frequency square. All these three methods should yield the same Fourier transformed accelerations(within numerical accuracy).

Source: http://prosig.com/wp-content/uploads/pdf/blogArticles/OmegaArithmetic.pdf

However, as you can see in the attached plot the magnitude of Fourier transformed acceleration is off by orders of magnitude. Can you please point out what I am doing wrong? I understand that the peak frequency is the same but the magnitudes and overall shape of plots is worrying me.

Here is my code

N =2**8
a =1
t=np.arange(N)
t = np.linspace(-100,100,N)
fs= 2
freq_normalized = rfftfreq(N)*fs

y = np.sin(a*t)
vel = np.cos(a*t)
acc = -np.sin(a*t)

fft_y = np.abs(rfft(y))
acc_y = [a*b**2 for a,b in zip(fft_y,freq_normalized)]

fft_v = np.abs(rfft(vel))
acc_v = [a*b for a,b in zip(fft_v,freq_normalized)]

acc_a = np.abs(rfft(acc))

acc_ifft_v = irfft(acc_v)
acc_ifft_y = irfft(acc_y)
acc_ifft = irfft(acc_a)

fig, ax = plt.subplots(figsize=(12, 8))
ax.spines['left'].set_position('center')
ax.spines['right'].set_color('none')
ax.spines['bottom'].set_position('center')
ax.spines['top'].set_color('none')
ax.spines['left'].set_smart_bounds(True)
ax.spines['bottom'].set_smart_bounds(True)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
plt.plot(t,y)
plt.ylabel(r'$ y $', fontsize=16)
plt.xlabel(r'$ t $', fontsize=16)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.show()

fig, ax = plt.subplots(figsize=(12, 8))
plt.title('Fourier spectrum from position')
plt.plot(freq_normalized,acc_y)
plt.ylabel(r'$ a_{\omega} $', fontsize=16)
plt.xlabel(r'$ \tilde \omega $', fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()

fig, ax = plt.subplots(figsize=(12, 8))
plt.title('Fourier spectrum from velocity')
plt.plot(freq_normalized,acc_v)
plt.ylabel(r'$ a_{\omega} $', fontsize=16)
plt.xlabel(r'$ \tilde \omega $', fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()

fig, ax = plt.subplots(figsize=(12, 8))
plt.title('Fourier spectrum from acceleration')
plt.plot(freq_normalized,acc_a)
plt.ylabel(r'$ a_{\omega} $', fontsize=16)
plt.xlabel(r'$ \tilde \omega $', fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()

The problem is that you are analyzing sine and cosine without using a smooth window e.g a Hanning window. Your phase range is from -100 to 100 giving 100 / $\pi$ = 31.83 periods.This means that you have a discontinuity at the end of the FFT time record. So you are not analyzing a sine or cosine as you believe. Using a Hanning window will suppress the discontinuities and give a more meaningful result.


Signal corrupted using peak (biquad) filter at low frequencies I'm trying to process an audio stream through a series of peak filters at static frequencies that range from 0 to nyquist. If I start the range from anything below about 25Hz, the audio gets corrupted in some way (I haven't looked at the values yet, so I don't know what the corruption looks like) that results in no sound. But making the range ~25Hz to nyquist is fine. This happens even when all the filters are at zero gain.

Is anyone aware of theoretical problems with using a low frequency peak filter? Maybe it's a problem when its bandwidth might cross over 0Hz? (Do standard biquad equations address that possibility?)

Looks a like a typical stability problem with the poles too close to the unit circle. A third octave wide parametric EQ filter at 25Hz sampled at 44.1kHz has a pole radius of > 0.999 so it's getting VERY close to the unit circle. The lower the frequency, the large the pole radius.The pole section alone has a peak gain for more than 110 dB !!!

With the poles so close to being unstable, you may run into all sorts of problem: numerical noise, noise amplification, coefficient rounding or truncation errors, filter topology (do NOT use Direct Form II or transposed Form I), etc.


Relation between bandwidth and data rate I have trouble understanding the following from wikipedia - Bandwidth (Signal Processing)

  A key characteristic of bandwidth is that any band of a given width can carry the same amount of information, regardless of where that band is located in the frequency spectrum

The higher the frequency the more the cycles and therefore shouldn't there be more information passing in a second?

A pure never changing sinusoidal carrier of a given frequency carries no information.  Changes or modulations to that carrier over time (as in: per second) can carry information.  The maximum rate of information is related to the bandwidth and S/N ratio, not the number of unmodulated carrier cycles.

Of course, it may be easier to implement certain modulation schemes with a higher carrier to bandwidth ratio.


Audio/speech detection, verification I would like to write a program to compare two '.wav' files, one of which is the original played file with human speech and the other is a recorded '.wav' file. The program I wish to write in python should be able to do the following:

Read the '.wav' files, compare them and say if speech exists or not. If yes, what is the match % and gain or loss if any.
Detect if there's any noise that's present in the recorded '.wav' file.
Calculate the delay between the original played file and the recorded file (playing and recording done using threads, so t=0 is the same for either files).

As far as requirement 3 is concerned, I was fairly successful in calculating the delay using correlation function that's offered by numpy as:

corr_files = numpy.correlate(orig_pitch, rec_pitch, 'full')
delay = int(len(corr_files) / 2) - numpy.argmax(corr_files)

For requirements 1 and 2, I happened to play around with an RMS based approach, but it didn't go well. What is the best way to achieve all the above requirements? Which algorithm suits the above requirements?

wav files generally have uncompressed data. That means you can get the exact waveform of the speech signal in digital form. You can tell whether the signal has some speech or not by calculating energy and using threshold to make the decision.

Again, the above technique may not work if there is noise along with speech.

Speech signal is quasi stationary signal and same speech will have different waveforms based on who is speaking, his pitch, his slang, his emotion and many more factors. So, it is very difficult to match two speech files.

The two problems that you wanted to solve

% of match between two files.
Detect noise.

are special fields in signal processing by itself. First one can be achieved by doing speech recognition on both the files and get the exact text for both files and then decide some rules to get % of match. For example, even though waveform of both speech are different but content is same, only with speech recognition you can 100% match.

coming to noise in the speech signal, there is again separate field in signal processing called as noise reduction. even though the process of noise reduction is used to reduce noise, one of the steps in noise reduction is to detect noise. You can use that step to detect noise.

But, both these topics are pure signal processing topics and i don't think you can achieve these with some simple function call or some simple python code.


Correct way to simulate Rayleigh fading for indoor environments I am trying to simulate indoor communications with all components included. But something disturbs me when it comes to Rayleigh fading part.

I am adding my rayleigh fading to the transmitted signal with that code:

rayleigh_h = 1/sqrt(2)*[randn(1, length(s_l) / OSR) + 1j * randn(1, length(s_l) / OSR)]; %Constructing rayleigh channel complex coefficients
rayleigh_h = repelem(rayleigh_h, OSR); %Coefficient repeating to apply rayleigh fading to all samples
s_l = s_l .* rayleigh_h;

s_l: the complex modulated signal

OSR: Over sampling ratio

If I need to clarify the code a bit: Every bit is being sampled by a number of samples named "Over sampling ratio". Since I want to add rayleigh fading symbol-by-symbol instead of sample-by-sample, all samples that belong to a symbol are multiplied by the rayleigh coefficient corresponding each bit.

Sorry for if the code clarification is a bit confusing.

Here is my real question.

The method I am using above obviously changing the amplitudes of the signal with multiplying the with complex gaussian coefficients at instants. However, at indoor locations, it is fairly correct to assume that mobility(speed) of the devices should be limited with the walking speed. Therefore changes with the signal, and the fading coefficients should change slowly, not that fast. I am trying to figure out how to simulate it.

Any help is appreciated. Please also advise me if I have misunderstood the concept of the rayleigh fading.

A discrete-time baseband model of a multi-propagation channel can be written as

$$y[m] = \sum_l h_l[m] x[m-l] + w[m]$$

where $l$ is the index of channel taps if the channel is modeled as a FIR $\{h_l, 0\leq l < L-1\}$. $w[m]$ is AWGN sample.

The channel tap $h_l[m] = \sum_i \alpha_i(m/W) \times \mathrm{sinc}(l-W\times\tau_i(m/W))$ where $i$ is the index of physical propagation path; $\alpha_i$ and $\tau_i$ are the gain and the delay of path $i$, respectively; $W$ is the bandwidth of signal $x(t)$.

The philosophy of this modeling is that due to the decay of $\mathrm{sinc}()$ function, $h_l[m]$ has significant contribution of the path $i$ having delay $\tau_i \in [l/W-1/2W,l/W+1/2W]$. Next the reason of Rayleigh modeling is that these number of that kind of path $i$ for $h_l$ is quite high and by the central limit theorem, $h_l$ is complex Gaussian.

The Rayleigh fading assumption is reasonable for the indoor environment because of the high number of scattering object.

The fact that you have used element-wise multiplication between s_l and rayleigh_h implies that your channel is one tap, i.e. $L=1$. One reason is that bandwidth $W \ll \textrm{delay spread } \tau_m$. Check it.

To model the fact that channel varies in time, you need to consider the notion of coherence time. The channel model you are trying to use is called Rayleigh Block fading. For a quick and dirty approach, just calculate the coherence time $T_c$ (based on the carrier frequency and the velocity), then divide it by the symbol time $1/W$ for the number of symbol that channel is assumed constant. I mean something like

nrep = OSR * Tc * W;
rayleigh_h = 1/sqrt(2)*[randn(1, length(s_l) / nrep) + 1i * randn(1, length(s_l) / nrep)];
rayleigh_h = repelem(rayleigh_h, nrep);
s_l = s_l .* rayleigh_h;

For more sophisticated model, I suggest using the rayleighchan object.


Amplitude calculation I've got an array of different values representing the y-values in a plot:

y-values = [0.07, 0.09, 0.04, 0.15, 0.14, 0.05, 0.02, 0.2, 0.18,...  ...0.06, 0.2, 0.17, 0.19, 0.21, 0.3]

I want to split it in 3 intervals:

the first one should return the values to the first big amplitude (here: [0.07, 0.09, 0.04, 0.15])
the second one should return the values to the last big amplitude (here: [0.14, 0.05, 0.02, 0.2, 0.18,...  ...0.06, 0.2])
the last one should contain the last values (here: [0.17, 0.19, 0.21, 0.3])

I have no idea how to calculate these threshold-amplitudes. Can anyone help me with this problem?

I have to implement this. I originally posted it on StackOverflow, but this is too mathematical as a programming subject, hope I can get some help here.

You have to be specific about your selection criteria. All you're "big amplitude" are local maxima but there are also other local maxima that don't qualify (why?). It looks like you want a local maximum that's above a certain threshold. Or maybe you only want the top 3?
Anyway, you can find local maxima by either differentiating and finding the zeros or my a simple search by looking at the neighbors of each sample. Then this can be combined with a threshold detector or a sorting algorithm.


Matching 2 audio files I'm doing research on how we can do matching between 2 audio files using Android to see if they have the same content or not. After doing lots of research, I have concluded that I must do 2 main things:

FFT
Correlation

Am I on the right track?

It really depends on which application you have...the algorithms for speech recognition are different from speech indentification for example. Based on what you need, just look for an open source library for Java. Getting into the theory is definitely desired but as a first step you can just look for a software block that has inputs/outputs, later you can drill into it.

Here is some info on speech recog: https://stackoverflow.com/questions/3018418/open-source-simple-speech-recognition-in-c-in-windows


How to adapt 4th order IIR filter coefficients from matlab to cmsis dsp biquad struct I'm designing an IIR 4th order filter to be implemented using CMSIS DSP on a cortex-M4 using a direct form II trannsposed structure. I get this coefficient from matlab's "butter" function:(consider Fs = 48kHz, bandpass Fc1 = 445 Hz, Fc2 = 561Hz)

b[5.70278412657461e-05   0   -0.000114055682531492   0   5.70278412657461e-05]
a[1   -3.97002047467733   5.9189048757951   -3.92762130007361   0.97875498919531]

I know i have to use two stages of a biquad cascade IIR filter and the according structure should receive coefficients in the follow way:

float32_t[10] = {b10, b11, b12, a11, a12, b20, b21, b22, a21, a22};

Question: Can some tell me how to translate matlab answer into CMSIS struct format?.

% design the filter as poles, zeroes and gain
[z,p,k] = butter(2,[445 561]*2/48000);
% turn into biquads
sos = zp2sos(z,p,k);
% eliminate a0
coeffs = sos(:,[1 2 3 5 6]);
% make a linear array
coeffs = [coeffs(1,:) coeffs(2,:)];


How to do the highest resolution FFT? I wish to do the FFT of streaming audio data coming into a computer. I am just interested in frequencies from 80 Hz up to 1000 Hz but with a frequency resolution of 3 to 4 Hz and no more.

Basically my question is , as my sampling frequency (Fs) is 44100, will my FFT window size have to be no less than 11025 (N) ?  (because spacing of FFT bins = Fs/N). Does this mean that I am completely limited to taking windows of data any higher than 4 times a second with my desired resolution.

Unfortunately I need to be able to capture windows of data 8 times a second and also have a resolution of 3 to 4 Hz.

Is this actually possible? Many thanks in advance.

The FFT size can be longer than the window size.  For instance you can use audio sample windows of 5512 (about 8 per second) and zero-pad them to length 16384 for each FFT.  If noise level is low enough and the separation between frequency peaks is great enough, the interpolation provided by this zero-padding can provide better frequency estimation and plot resolution.

The "highest resolution" FFT reasonable depends on the signal-to-noise ratio as well as the window length, as after a certain point, interpolation and frequency estimation methods will be affected as much by nearby spectral noise as by the spectral peaks that one is trying to estimate at "high resolution".


Measuring the similarity of two waveforms I have two approximately sine truncated waveforms of different amplitudes and I want a numerical measure of how similar they are to each other. I am sending a sonic pulse train whose multiple reflections I am receiving. I want to know which of those reflections most closely correlate with each other in terms of wave shape.
Any standard methods?

It seems the term (and process) I am looking for is "autocorrelation".


Gaussian filter looks pixelated What does it mean when an image looks pixelated after a Gaussian filter (Multivariate Gaussian Distribution filter) is applied (in context of mu, Sigma, and the meshgrid created), and what parameters can one alter to minimize this effect and make the image look smoother?

Sample code (MATLAB):

d1 = 20
d2 = 20
d3 = 20
mu = [0 0 0];
Sigma = 30*[1 0 0; 0 1 0; 0 0 1];
[X1,X2,X3] = meshgrid(linspace(-10,10,d1)', linspace(-10,10,d2)', linspace(-10,10,d3)');
X4 = [X1(:) X2(:) X3(:)];
p = mvnpdf(X4, mu, Sigma);
F = reshape(p,d1,d2,d3);

I am convolving the filter above with a 3-D matrix, a roughly 90*60*100 matrix.

An image shouldn't look pixelated after filtered by a Gaussian (lowpass) kernel... instead look at the following figure:


Energy of a complex exponential I want to obtain the energy of this signal $$x_i(t)= e^{j2\pi(2i-1)f_0t}$$ where $f_0 = 1$ Hz, $i= 1,...,4$ and $j$ is a complex. I think I need to use the Fourier transform or the Parseval theorem but I'm not sure how to apply it.

As mentioned in the comments you can just take the integral of the squared magnitude.
Notice that the magnitude of $|e^{a+bj}| = |e^a||e^{bj}| = e^a$ for $a,b \in \mathbb{R}$, comparing with your function we see that $|x_i(t)|=1$. We conclude that $x_i(t)$ is a power signal, not an energy signal, i.e. you must define a time interval to integrate. Then the energy will be $\int_{t_0}^{t_f} |x_i(t)|^2  \mathrm{d}t = (t_f - t_0)$. I hope this helps.


Does anyone know how to overlay a equalized CDF plot with equalized histogram in matlab? I am willing to overlay a equalized cdf with equalized histogram in matlab.
But it's really hard to find any recommends.
so if you can, would you please let me know about that?

You can use 'hold on' with a histogram:

[n,x]=hist(rand(1,100),20);  % find the histogram
bar(x,n./sum(n))             % plot normalized histogram
hold on                      % don't delete current plot
plot([0 1],[0.1 0.1])        % plot a uniform pdf
axis([0 1 0 0.12])           % make room around the plot

Another example using a normal distribution, which shows how to normalize the histogram so that its area is one:

r=randn(1,1000);      % normal random variable
[n,x]=hist(r,25);
delta=x(2)-x(1);
n=n./(delta*sum(n));  % n is normalized so that its area is 1
bar(x,n)
hold on
d=-4:0.1:4;
plot(d,(1/sqrt(2*pi))*exp(-(d.*d)./2),'r')  % plot of the normal pdf


Noise Characterization i am working on characterization of noise in an experimental signal from an accelerometer, mounted on a beam. Normally, white gaussian noise is considered by several researchers to add noise in the simulated model. I want to find out, what is the noise type in real scenarios, whether it is gaussian or non-gaussian.

There are 2 properties which are important to know about your noise:

Its Self Correlation
What usually called Auto Correlation. This will tell you the relation between 2 noise samples which are apart (On time) from each other (I'm assuming the signal is stationary, namely its properties doesn't change over time).
Best thing to hope for is White Noise.
Its Distribution
Sometimes, in order to be able to do estimations of the data one should also know the distribution of the noise. The most wide spread model is Gaussian Noise as you wrote. Yet on Image Processing it could be Poisson Noise.

If you have access to "Pure" noise samples analyzing those properties of the noise is really easy.
For the first you just calculate the Auto Correlation Function.
For testing Gaussianity / Normality Test (Reference - Tests of Gaussianity) things gets a little tricky.
As you can see in the Wikipedia Page of Normality there are many tests.
Though poor man's (Yet effective) test would be just looking at the Empiric Histogram and compare it to the model (Curve of Normal PDF based on the data Mean and Variance).


Consider the integrator and check for time invariance I'm trying to check the following system for for time-invariance.

$$ y(t) = \int_0^t x(\lambda) d\lambda$$

Please explain me why it is time-invariant.

Note that there are different kinds of (causal) integrators. The one in your question with a fixed finite lower integration limit is not time-invariant, as you've shown yourself (in the comments).

Other variants are

$$y(t)=\int_{-\infty}^tx(\tau)d\tau\tag{1}$$

and

$$y(t)=\int_{t-T}^tx(\tau)d\tau\tag{2}$$

Both of these integrators are time-invariant, which is straightforward to show.


Relationship between interpolation accuracy, impulse response and frequency response Global interpolation or sinc interpolation is an ideal filter since its frequency response is a rect function.
The impulse response of this filter is the sinc function (same as the coefficients of the interpolator).

Local interpolation methods are constructed by applying some window to the sinc function.
For example, according to https://ccrma.stanford.edu/~jos/Interpolation/Relation_Lagrange_Interpolation_Windowed.html, for uniformly spaced samples and finite $N$ , Lagrange interpolaton is equivalent to windowed sinc interpolation using a binomial window.
The impulse response of this filter are the coefficients of the interpolator.

I am interested to know what is the relationship between the impulse response, accuracy of the interpolation and frequency response of the filter/local interpolator, especifically for Lagrange interpolation.
I cannot see this clearly.

any interpolation polynomial that goes through the two adjacent samples to the left and right of the area being interpolated will result in an effective impulse response that is a piecewise-polynomial function going through zero at all integers except zero (like the sinc function does).

e.g. linear interpolation:

third-order lagrange:

third-order hermite:

note all of these effective impulse responses go through the $(x,y)$ points of $(0,1)$ and $(n,0)$ for any other integer $n \ne 0$, just like the sinc function does.

so you can represent that interpolation function as something times a sinc function and the "something" can be obtained by dividing that interpolation function, point-by-point, with the sinc function.  the only division by zero is when the numerator is also zero and that can be, what we like to call in the hand-waving department, a "removeable singularity".  (maybe we can find a decent limit at the $\tfrac00$ divisions.)

any interpolation that forces the smooth interpolated function to go through the discrete points can always be described as a windowed-sinc function.  all you have to do is divide to get the effective window.


Software Mixing of Multiple Audio Streams From Network TLDR: How do I mix network streams when streams might be unreliable?
I am working on a project which will allow musicians to practice/perform together over the internet with low latency. (Think JackTrip but more user friendly).
I am receiving raw PCM audio data via UDP over the network. As these come in, I identify who the sender is and pipe that stream into that particular person's buffer. I have a separate buffer for each person participating.
Periodically I then take one "chunk" at a time from each buffer, mix them together and output to the speakers of the device.
The issue I'm having is knowing when to take that "chunk". Some clients go offline and their buffers run empty. I cannot just wait for their buffer to refill as this will slow down the other clients.
Thoughts I had were to look at all the buffers at a particular point in time and only process a chunk when at least half of the buffers are populated.
Any ideas?
Other details:

My buffers are 128 samples big.
I am using a sample rate of 44100
I am using 16 bit signed integer pcm
Streams are MONO
I am using PortAudio and C++.
I am using ASIO audio drivers.

Generally you will need some sort time sync protocol.
Typically you want create locally one circular buffer for each network channel. As soon as a packet comes in, you shove it into circular buffer ort that channel at the write positions.
Playback is clocked by your local D/A clock. As soon as the D/A buffer is running empty, you grab a frame from all the client circular buffers at the current read position, mix them together and put them into the D/A buffer.
There are two things you need to manage.

Drop out: That happens if you need a new playback frame and the circular buffer of a client is empty OR a packet got dropped (they should be numbered to detect this). There isn't a whole lot you can do about. You can play zeros or try to extrapolate from the previous buffer A drop out means that your network can't keep up.
Clock drift. Even if everyone is playing at nominally 44.1kHz, no two clocks are identical and over time they will drift apart. To fix these you need to a) measure the clock drift and b) correct it. a) can be done through either a global time reference (something like NTP), by averaging over packet arrival times, or monitoring buffer levels. b) will require some sort of sample rate conversion: the most crude option would be to occasionally drop a frame insert a zero frame.

This is age old trade off: the lower the latency, the higher your drop out rate will be. So dialing this in "optimally" is the key here. It's easy enough to adjust: simply make your circular buffers smaller/bigger. However finding the optimum is hard since it depends very much on current network conditions so in practice this really should be adaptive.


Expression for Frequency-shifted Gaussian pulse I am doing some research on UWB radars which transmit frequency-shifted Gaussian pulses. These pulses are given by:
$$
g(t)  = p(t)cos(2\pi f_{c}t) = V_{TX}exp(-\frac{t^{2}}{2\tau^2})cos(2\pi f_{c}t)
$$
where $f_{c}$ is the center or carrier frequency, $V_{TX}$ is the peak amplitude and $\tau$ determines the Bandwidth of the pulse (you can think of $\tau$ as the standard deviation of $p(t)$. This expression in some papers is equivalent to a higher order derivative of the Gaussian function $p(t)$, such that:
$$
g(t) = \frac{d^{n}}{dt^{n}} \left( V_{TX }exp(-\frac{t^{2}}{2\tau^2}) \right)
$$
where $n$ denotes the derivative order.
Does anyone know how these two expressions are related? Is there a way of obtaining one from the other?

Well the expression for the $\frac {d^n}{dt^n} p(t) $ is $V_n(t) p(t)$, where $V_n(t)$ is a polynomial of degree $n$, and $V_n(t)$ may be determined recursivelyy
$$  V_{n+1}(t)e^{-t^2/(2\tau^2)} = \frac{d}{dt}\left(V_{n}(t)e^{-t^2/(2\tau^2)}\right) $$
Solving for $V_{n+2}(t)$ we have
$$ V_{n+1} = \frac{d}{dt} V_{n}(t) - \frac{t}{2\tau^2} V_n(t)$$
In order to link the two forms of the equation we need
$$\begin{eqnarray}
 g(t) &=&  p(t)cos(2\pi f_{c}t) \\
      &=& V_{TX}exp(-\frac{t^{2}}{2\tau^2})cos(2\pi f_{c}t) \\
      &=& \frac{d^{n}}{dt^{n}} \left( V_{TX }exp(-\frac{t^{2}}{2\tau^2}) \right) \\
      &=& p(t)V_n(t)
\end{eqnarray}$$
Thus $cos(2\pi f_c t) = V_n(t)$, the right hand side is a polynomial. As $n \to \infty$, the are the Hermite polynomials have some asymptotic expansions in terms of cossines. So the two may assume approximate values for certain values of $f_c$ and $\tau$.
Numerically evaluating the derivaties we see some resemblance with gausiand windowed cosine.
function hermite;
  tau=pi;
  t = tau * linspace(-5, 5, 1000);
  for order = 4:2:20
    subplot(3,3,order/2-1)
    a = -1/(2*tau^2);
    P = V_coefs(order, a);
    plot(t, polyval(P(end:-1:1), t).*exp(a*t.^2));
    title(num2str(order, 'n=%d'))
  end

end

function coefs = V_coefs(order, a)
  coefs = 1;
  for ii = 1:order
    coefs = [(1:(ii-1)).*coefs(2:end), 0, 0] + [0, coefs .* (2 * a)];
  end
end


[Python ]How can I improve my 1D Kalman Filter estimate? I have written the following code to smooth an (almost) linear function:

import numpy as np
import matplotlib.pyplot as plt

# intial parameters
n = 500
sz = (n) # size of array
t = np.linspace(0,10,n)
x = t + 2*t**2 # truth function
z = x + np.random.normal(0,7.5,size=sz) # noise

Q = 1e-4 # process variance

# allocate space for arrays
xhat=np.zeros(sz)      # a posteri estimate of x
P=np.zeros(sz)         # a posteri error estimate
xhatminus=np.zeros(sz) # a priori estimate of x
Pminus=np.zeros(sz)    # a priori error estimate
K=np.zeros(sz)         # Kalman gain

R = 0.1**2 # estimate of measurement variance

xhat[0] = 0.0
P[0] = 1.0

for i in range(1,n):

  xhatminus[i] = xhat[i-1]
  Pminus[i] = P[i-1]+Q

  K[i] = Pminus[i]/( Pminus[i]+R )
  xhat[i] = xhatminus[i]+K[i]*(z[i]-xhatminus[i])
  P[i] = (1-K[i])*Pminus[i]

plt.figure()
plt.plot(t, z,'k+',label='noisy measurements')
plt.plot(t, xhat,'b-',label='a posteri estimate')
plt.plot(t,x,color='g',label='truth value')
plt.legend()
plt.xlabel('t')
plt.ylabel('f(t)')

plt.show()

I have chosen x[0] = 0 because I know that's a good approx and P[0] = 1 because I didn't know better.

Plot of the estimate, noise and the true function:

I know that Kalman filters are not meant for non-linear functions like mine!

I have played around with both Q and R to try to find a better estimation of the true function.My function converges very slowly.  Is there a way to find a "better" Q and R?

I believe ZHUANG's answer is correct, so I went searching for some source code in Python that implemented it.  The closest I could find was a 2D example that uses velocity as well.  The reference is at the top of the listing.  I have reformatted and restructured the code to make it more readable to me and likely more efficient.  I also adapted it to the OPs function.

I am a newbie to Kalman filters.  I am using this as my learning tool.  It can probably be improved upon greatly.  I am posting it hoping that it will be helpful to the OP as well as others.

Ced

# Derived from:
# https://stackoverflow.com/questions/13901997/kalman-2d-filter-in-python
#   and
# https://dsp.stackexchange.com/questions/48343

import numpy as np
import matplotlib.pyplot as plt

#=============================================================================
def main():

#---- Set Parameters

        N = 50

        x = np.matrix( '0. 0. 0. 0.' ).T
        P = np.matrix( np.eye( 4 ) ) * .0001000 # initial uncertainty

        R = 0.1**2

#---- Set up Helper Matrices

        F = np.matrix( '''
            1. 0. 1. 0.;
            0. 1. 0. 1.;
            0. 0. 1. 0.;
            0. 0. 0. 1.
            ''' )

        H = np.matrix( '''
            1. 0. 0. 0.;
            0. 1. 0. 0.''' )

        motion = np.matrix( '0. 0. 0. 0.' ).T

        Q = np.matrix( np.eye( 4 ) )

        m = np.matrix( '0. 0.' ).T

#---- Create the Test Case

#        true_x = np.linspace( 0.0, 10.0, N )
#        true_y = true_x**2

#        observed_x = true_x + 0.05 * np.random.random( N ) * true_x
#        observed_y = true_y + 0.05 * np.random.random( N ) * true_y

        sz = (N) # size of array

        true_x = np.linspace( 0.0, 10.0, N )
        true_y = true_x + 2 * true_x**2

        observed_x = true_x
        observed_y = true_y + np.random.normal( 0, 7.5, size=sz )

#---- Run the Sequence

        kalman_x = np.zeros( N )
        kalman_y = np.zeros( N )

        for n in range( N  ):
            m[0] = observed_x[n]
            m[1] = observed_y[n]
            x, P = kalman( x, P, m, R, motion, Q, F, H )
            kalman_x[n] = x[0]
            kalman_y[n] = x[1]

#---- Plot the Results

#        plt.plot( observed_x, observed_y, 'ro' )
#        plt.plot( kalman_x, kalman_y, 'g-' )
#        plt.show()

        plt.figure()
        plt.plot( true_x, observed_y,'k+',label='noisy measurements')
        plt.plot( true_x, kalman_y,'b-',label='a posteri estimate')
        plt.plot( true_x, true_y,color='g',label='truth value')
        plt.legend()
        plt.xlabel('t')
        plt.ylabel('f(t)')

        plt.show()

#=============================================================================
#def kalman( x, P, measurement, R, motion, Q, F, H ):
def kalman( x, P, m, R, motion, Q, F, H ):
        '''
        Parameters:
        x: initial state
        P: initial uncertainty convariance matrix
        measurement: observed position ( same shape as H * x )
        R: measurement noise ( same shape as H )
        motion: external motion added to state vector x
        Q: motion noise ( same shape as P )
        F: next state function: x_prime = F * x
        H: measurement function: position = H * x

        Return: the updated and predicted new values for ( x, P )

        See also http://en.wikipedia.org/wiki/Kalman_filter

        This version of kalman can be applied to many different situations by
        appropriately defining F and H
        '''

#---- UPDATE x, P based on measurement m
#            distance between measured and current position-belief

        y = m - H * x

        S = H * P * H.T + R  # residual convariance
        K = P * H.T * S.I    # Kalman gain

        x = x + K * y
        I = np.matrix( np.eye( F.shape[0] ) ) # identity matrix
        P = ( I - K * H ) * P

#---- PREDICT x, P based on motion

        x = F * x + motion
        P = F * P * F.T + Q

#---- Exit

        return x, P

#=============================================================================
main()


Dual-Tone Multi Frequency signalling Does anybody know why we need 2 different tone signals: one high frequency and one low frequency to create a tone burst in Dual-Tone Multi Frequency signalling?

For example, to create a tone burst of 3, we need a high frequency $1477 \mathtt{Hz}$ and a low frequency $697 \mathtt{Hz}$. Does anyone know the reason for that?

Essentially, the reason why you need two tones is to ensure that a normal human voice can't replicate the tones.  As such, by simultaneously generating a tone from a high-frequency group and low-frequency group, it is highly improbable that a human voice can replicate a sound from such differing ends of the frequency spectrum at the same time.

If we were to adopt single tone dialing, wouldn't it be annoying if you were to talk to someone, your voice matches one of the dialing frequencies and it would randomly dial that number while talking?  It would certainly annoy me!

The Dual-Tone Multi Frequency tones (or DTMF for short) for a particular key is split up in the following grid:

1     |    2     |   3    |    A    |   697Hz
4     |    5     |   6    |    B    |   770Hz
7     |    8     |   9    |    C    |   852Hz
*     |    0     |   #    |    D    |   941Hz
--------------------------------------
1209Hz   1336Hz    1477Hz    1633Hz

For example, by pushing 1, you would generate two tones: 697 Hz and 1209 Hz from the low and high frequency ranges respectively.  In your post, you have already identified the frequencies used when dialing the number 3.  A,B,C,D are probably something you haven't seen before.  A,B,C,D were actually used as keys for the US Military's Autovon phone system.  They were used for:

A - Flash Override
B - Flash
C - Immediate
D - Priority

Pushing one of these keys gave your call priority with Flash Override being the highest.  The Autovon system has been disbanded for a long time, but many non-telephone applications still use those keys today.

One question that you may ask is how they decided to come up with these frequencies?  The frequencies were chosen to avoid harmonics - as in no frequency being a multiple of another.  Also, the difference between any two frequencies doesn't equal any of the other frequencies and the sum of any two frequencies isn't equal to any of the frequencies too.  The frequencies were generated within a +/- 1.5% tolerance of each other, the highest frequency must be as loud as the lowest frequency and as much as 4 dB louder.  This level difference is referred to as a "twist".


Generate the Matrix Form of 2D Convolution Kernel We know that a convolution can be replaced by a multiplication with a Toeplitz / Circulant Matrix. Meaning, assume I have convolution kernel $ h $ and matrix $ I $ (Of size $ m \times m $ for example), then there is a matrix $ H $ of size $ m^2 \times m^2 $ such that $ h \ast I $ is the same as $ H I^{cs} $ Where cs for column stacked image resized to $m \times m$ (reshape(H*I(:), [m, m])).

My question is - how to construct $H$? I know that every row of $H$ is responsible for single entry of the resulting multiplication (H(1,:)*I(:) gives the (1,1) entry of the convolution result), so I can try to work very hard to calculate each entry of $H$ (e.g. if the kernel $h$ is of size $3 \times 3$, than the convolution result in entry (i,j) is:

$$ h_{11}\cdot I(i-1,j-1) + h_{12}\cdot I(i-1,j) + h_{13}\cdot I(i-1,j+1) + h_{21}\cdot I(i,j-1) + h_{22}\cdot I(i,j) + h_{23}\cdot I(i,j+1) + h_{31}\cdot I(i+1,j-1) + h_{32}\cdot I(i+1,j) + h_{33}\cdot I(i+1,j+1) $$

except for near the edges of $I$. so it's somewhat possible with a lot of hard work to construct $H$ such that when multiplying by $I^{cs}$ I will get that result). But is there a 'generic' way for this construction? Or perhaps a made MATLAB function?

I am aware of the convmtx2() function in MATLAB, but the resulting matrix is not in the proper dimensions (not even square).

I created a function to create a Matrix for Image Filtering (Similar ideas to MATLAB's imfilter()):

function [ mK ] = CreateImageFilterMtx( mH, numRows, numCols, operationMode, boundaryMode )
% ----------------------------------------------------------------------------------------------- %
% [ mK ] = CreateImageFilterMtx( mH, numRows, numCols, operationMode, boundaryMode )
% Generates an Image Filtering Matrix for the 2D Kernel (The Matrix mH)
% with support for different operations modes (Convolution / Correlation)
% and boundary conditions (Zeros, Symmetric, Replicate, Circular). The
% function should match the use of MATLAB's 'imfilter()' with the same
% parameters.
% Input:
%   - mH                -   Input 2D Convolution Kernel.
%                           Assumed to have odd dimensions.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numRows           -   Number of Rows.
%                           Number of rows in the output convolution
%                           matrix.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - numCols           -   Number of Columns.
%                           Number of columns in the output convolution
%                           matrix.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - operationMode     -   Operation Mode.
%                           Sets whether to use Convolution or Correlation
%                           for the operation mode.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2}.
%   - boundaryMode      -   Boundary Condition Mode.
%                           Sets the boundary condition mode for the
%                           filtering. The options are - Zeros, Symmetric,
%                           Replicate and Circular.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, 4}.
% Output:
%   - mK                -   Convolution Matrix.
%                           The output filtering matrix. Multiplying in
%                           the column stack form on an image should be
%                           equivalent to applying 'imfilter()' on the
%                           image.
%                           Structure: Matrix (Sparse).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'convmtx2()' - https://www.mathworks.com/help/images/ref/convmtx2.html.
% Remarks:
%   1.  The height and width of 'mH' are assumed to be odd number. In case
%       either or both are even the user should pad the kernel with zeros
%       (Either a row, column or both). according to the anchor of the
%       kernel the user do the padding pre or post the kernel.
% TODO:
%   1.  Refactor the code to share the common operations of different
%       boundary modes.
% Release Notes:
%   -   1.0.001     30/12/2019  Royi Avital
%       *   Fixed some bugs related to using 'numCols' instead of 'numRows'
%           in the calculation of 'pixelShift' for the cases 'jj + ll >
%           numCols' and 'jj + ll < 1'.
%   -   1.0.000     16/01/2018  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

OPERATION_MODE_CONVOLUTION = 1;
OPERATION_MODE_CORRELATION = 2;

BOUNDARY_MODE_ZEROS         = 1;
BOUNDARY_MODE_SYMMETRIC     = 2;
BOUNDARY_MODE_REPLICATE     = 3;
BOUNDARY_MODE_CIRCULAR      = 4;

switch(operationMode)
    case(OPERATION_MODE_CONVOLUTION)
        mH = mH(end:-1:1, end:-1:1);
    case(OPERATION_MODE_CORRELATION)
        % mH = mH; %<! Default Code is correlation
end

switch(boundaryMode)
    case(BOUNDARY_MODE_ZEROS)
        mK = CreateConvMtxZeros(mH, numRows, numCols);
    case(BOUNDARY_MODE_SYMMETRIC)
        mK = CreateConvMtxSymmetric(mH, numRows, numCols);
    case(BOUNDARY_MODE_REPLICATE)
        mK = CreateConvMtxReplicate(mH, numRows, numCols);
    case(BOUNDARY_MODE_CIRCULAR)
        mK = CreateConvMtxCircular(mH, numRows, numCols);
end

end

function [ mK ] = CreateConvMtxZeros( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                % Pixel Index Shift such that pxIdx + pxShift is the linear
                % index of the pixel in the image
                pxShift = (ll * numRows) + kk;

                if((ii + kk <= numRows) && (ii + kk >= 1) && (jj + ll <= numCols) && (jj + ll >= 1))
                    vCols(elmntIdx) = pxIdx + pxShift;
                    vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);
                else
                    vCols(elmntIdx) = pxIdx;
                    vVals(elmntIdx) = 0; % See the accumulation property of 'sparse()'.
                end
            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);

end

function [ mK ] = CreateConvMtxSymmetric( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                % Pixel Index Shift such that pxIdx + pxShift is the linear
                % index of the pixel in the image
                pxShift = (ll * numRows) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - (2 * (ii + kk - numRows) - 1);
                end

                if(ii + kk < 1)
                    pxShift = pxShift + (2 * (1 -(ii + kk)) - 1);
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - ((2 * (jj + ll - numCols) - 1) * numRows);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + ((2 * (1 - (jj + ll)) - 1) * numRows);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);

end

function [ mK ] = CreateConvMtxReplicate( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                % Pixel Index Shift such that pxIdx + pxShift is the linear
                % index of the pixel in the image
                pxShift = (ll * numRows) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - (ii + kk - numRows);
                end

                if(ii + kk < 1)
                    pxShift = pxShift + (1 - (ii + kk));
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - ((jj + ll - numCols) * numRows);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + ((1 - (jj + ll)) * numRows);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);

end

function [ mK ] = CreateConvMtxCircular( mH, numRows, numCols )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

numElementsImage    = numRows * numCols;
numRowsKernel       = size(mH, 1);
numColsKernel       = size(mH, 2);
numElementsKernel   = numRowsKernel * numColsKernel;

vRows = reshape(repmat(1:numElementsImage, numElementsKernel, 1), numElementsImage * numElementsKernel, 1);
vCols = zeros(numElementsImage * numElementsKernel, 1);
vVals = zeros(numElementsImage * numElementsKernel, 1);

kernelRadiusV = floor(numRowsKernel / 2);
kernelRadiusH = floor(numColsKernel / 2);

pxIdx       = 0;
elmntIdx    = 0;

for jj = 1:numCols
    for ii = 1:numRows
        pxIdx = pxIdx + 1;
        for ll = -kernelRadiusH:kernelRadiusH
            for kk = -kernelRadiusV:kernelRadiusV
                elmntIdx = elmntIdx + 1;

                % Pixel Index Shift such that pxIdx + pxShift is the linear
                % index of the pixel in the image
                pxShift = (ll * numRows) + kk;

                if(ii + kk > numRows)
                    pxShift = pxShift - numRows;
                end

                if(ii + kk < 1)
                    pxShift = pxShift + numRows;
                end

                if(jj + ll > numCols)
                    pxShift = pxShift - (numCols * numRows);
                end

                if(jj + ll < 1)
                    pxShift = pxShift + (numCols * numRows);
                end

                vCols(elmntIdx) = pxIdx + pxShift;
                vVals(elmntIdx) = mH(kk + kernelRadiusV + 1, ll + kernelRadiusH + 1);

            end
        end
    end
end

mK = sparse(vRows, vCols, vVals, numElementsImage, numElementsImage);

end

The code was validated against MATLAB imfilter().

Full code is available in my StackOverflow Q2080835 GitHub Repository.


Effect of changing sample rate, window duration and zero padding on DTFT and DFT Let $T$ be the window duration, $N$ be the DFT size, $F_s$ be the sample rate, and $F_{max}$ be the frequency of the highest bin. In the context of image below:

halving the $F_s$ (keeping $T$ constant and N constant halving N) would remove the right half of the bins (halving the maximum frequency)
halving $T$ would change the shape of the DTFT, halving the number of lobes. Without zero padding this would halve $N$, but $F_{max}$ would remain constant and so every second bin would be removed.
doubling $N$ by zero padding the DFT would keep the DTFT and $F_{max}$ constant but would double the number of bins (i.e. interpolation of the DTFT).

Is my understanding sound?

Edit: No. I had overlooked aliasing when halving $F_s$

halving the $F_s$ (keeping $T$ and $N$ constant) would remove the right half of the bins (halving the maximum frequency)

No. Number of bins remain same because you have not changed $N$. What will change is the resolution. Each bin will represent $F_s/(2N)$ instead of $F_s/N$. Your DFT will not represent analog frequencies from $F_s/2$ to $F_s$ That is why you misunderstood as 'remove right half of the bins'.

  halving $T$ would change the shape of the DTFT, halving the number of lobes. Without zero padding this would halve $N$, but $F_{max}$ would remain constant and so every second bin would be removed

I am assuming that when you halve $T$, you are also reducing DFT to $N/2$. If you reduce FFT size with no zero padding, yes, you are dropping every second bin. Because the resolution now increases to $F_s/(N/2)$. For a period signal, with $N/2$ multiple of period, you will still see only 2 peaks for the FFT because the lobe width also has doubled, and the zero crossing appear exactly at multiples of $F_s/(N/2)=2 F_s/N$.

  doubling N by zero padding the DFT would keep the DTFT and Fmax constant but would double the number of bins

Yes, zero padding would increase the resolution of DFT. $F_s/(2N)$ is the resolution as compared to earlier $F_s/N$. But remember by zero padding you are not adding any new information about the signal. It just reveals more of the DFT you already have.


Straight lines with an edge detection algorithm I've been thinking of an algorithm for object recognition but it would relies heavily on straight non-noisy lines and as far as I know this is hard to obtain with edge detection algorithms. What's the best edge detection algorithm to obtain lines as straight as possible?

You're probably looking for the Hough transform or one of it's extensions.

The simplest version of this transform is linear and appropriate for detecting straight lines.

In the transformed space (Hough space), angles and distances are found as points where curves intersect.

Libraries for calculating the Hough transform exist in

C++ - OpenCV (Has Python bindings)
Matlab - hough in the image processing toolbox
Python - Scikits-image

Related Questions and Answers

See this previous answer for more help on understanding the Hough transform.


How compute frequency of signal from zero crossing? I need to compute frequency of a signal from zero crossing.It's possible to compute to frequency from zero crossing of a signal.
Plz anybody help me about this with matlab .

From your signal crossing threshold it is fairly straightforward to calculate its period. When you know the period then calculating of corresponding frequency is very simple (its inverse). This answer: enter here answers to similar question.


Given the 4 corners of two QR codes on a photographed sheet of paper, can we get a deskewed rectangular image? I made a photo of 2 sheets of paper with a phone camera, each of them cropped at the:

(x1, y1) top left corner of top QR code
(x2, y2) bottom right corner of bottom QR code

I hoped that the x, y coordinate of checkboxes would be nearly the same. In fact it's not, because of the angle of the camera, skewness, etc.

Question: How to have a better deskewing by using 8 corners:

the 4 corners of top QR code  (detected automatically by this method): (x1, y1), ..., (x4, y4)
the 4 corners of bottom QR code: (x5, y5), ..., (x8, y8)

?

Note: links of the 2 images: https://i.stack.imgur.com/SR8hz.jpg and https://i.stack.imgur.com/XsxrT.jpg

The solution is to use an "homography" / aka "projective transformation" (see this PDF, page 16).

Here is a working code showing how to do it with Python + OpenCV.