Quality of experience engineer I have noticed certain overlap between digital signal processing and quality assurance for jobs that require expertise in evaluating the quality of the transmitted signals on telecommunication systems like VoIP or IPTV. I have two examples:

https://www.usajobs.gov/GetJob/ViewDetails/511503700
https://highfive.com/careers/positions/?gh_jid=202317

For me, however, electronics engineer or quality assurance engineer are not terms that describe enough someone who is able to understand ITU-T standards (MOS, PSQM, PEVG, PESQ, POLQA, MUSHRA), quality of experience (QoE), network terminology (packet loss, bit rate, throughput, transmission delay, availability, jitter), DSP techniques (frequency analysis, FIR/IIR, sampling and quantization), and software programming.

Is there another name for this position? Is there a career path defined for someone who aims to work in this industry? Is this job profile something required nowadays by other companies?

PS: I wonder also if an expert like this can also participate in the creation of new algorithms and the evaluation of the quality of the architecture of DSP systems, as described here: https://www.eetimes.com/document.asp?doc_id=1275520

In my experience in the defense industry, the quality assurance engineers are mostly just “box checkers” that have limited understanding of how all of the individual facets of the algorithms and systems work. We have design engineers, integration engineers, and then the QA engineers. I have never seen an integration engineer or QA engineer get involved in any sort of design; the skill sets for each position are very, very different. It’s one thing to know what something is, but it’s something else entirely to be able to use it in a creative way to solve a problem.

Again, this is just my experience in the defense industry within the US. Other countries and industries may be a bit different though!


Identify the frequency spectrum of an audio file by using matlab I need to identify the first few harmonics of a audio file by cancelling the noise of the specific signal.
So I have partitioned the signal into 1024 Hz sample blocks and then took the FFT values of each block separately and finally add all FFT values and plot the waveform. Here is my matlab code. Somebody help me with the code since i cannot add the first block fft into total fft. Get an error. Please somebody help me

clc
clear
[signal,fs] =audioread('Female_55.wav'); %Only for plotting
t=0:1/fs:(length(signal)-1)/fs;          %Only for plotting

fs0=fs/1024                        %number of frequency blocks
t1=length(signal)/fs;              %length of the waveform

s=1;    %starting value
s1=1;   %ending value
c1=1;

wavefft_tot= zeros(1025,1);

for blocks = 1:fs0

[signal1,fs1] =audioread('Female_55.wav',[s*c1 s1*1024]);
wavefft_next=abs(fft(signal1));
wavefft_prv=wavefft_next;

n=length(signal1)-1;  %plotting information
f=0:fs1/n:fs1;        %plotting information

if(s==1) && (c1==1)
    c1=1;
    s=1024;
    s1=s1+1;
    wavefft_tot= wavefft_next;
else
   s1=s1+1;
   c1=c1+1;
   wavefft_tot=wavefft_tot+wavefft_prv;

    end

    end

    figure(1);
    plot(f,wavefft_tot);
    xlabel('Frequency in Hz');
    ylabel('Magnitude');
    title('The Wave FFT');

It looks like you are trying to recreate the functionality in spectrogram. Spectrogram will give you a short-time Fourier transform using whatever blocksize you specify.

Much of your code looks like it is re-inventing this.


How do you relate imaginary numbers with phase shift? How to imagine this? So I found this answer on one of the questions. I wanted to ask how the complex part can be related with the phase shift. I wanted to ask this on the same question but it didn't allow me. Could anyone please help me understand this concept. Even if I imagine this as a spiral with the value of the imaginary part changing, the starting time of the wave remains the same doesn't it? So how is it related to phase shift?

"Negative frequency doesn't make much sense for sinusoids, but the Fourier transform doesn't break up a signal into sinusoids, it breaks it up into complex exponentials:

$F(\omega)=\int_{-\infty}^{\infty} f(t)e^{−jωt}dt$

These are actually spirals, spinning around in the complex plane:

complex exponential showing time and real and imaginary axes

Spirals can be either left-handed or right-handed (rotating clockwise or counterclockwise), which is where the concept of negative frequency comes from. You can also think of it as the phase angle going forward or backward in time.

In the case of real signals, there are always two equal-amplitude complex exponentials, rotating in opposite directions, so that their real parts combine and imaginary parts cancel out, leaving only a real sinusoid as the result. This is why the spectrum of a sine wave always has 2 spikes, one positive frequency and one negative. Depending on the phase of the two spirals, they could cancel out leaving a purely real sine wave or a real cosine wave or a purely imaginary sine wave, etc."

Paint a red dot on the blue complex spiral at t=0, where the value is 1.0.  Now roll the spiral upside-down until the red dot is at -1.0, still at t=0.  This roll would be equivalent to a phase shift of Pi.  Note that the real part no longer corresponds to a cosine wave that starts at t=0.

But the real part of the rolled spiral does correspond to a cosine wave that starts (with value 1.0) at +Pi or -Pi, which is a different starting time, contradicting one of your assumptions above.

Thus phase is related to starting time.


Why can't the callee whom I am on call with hear the sound coming from my phone's speaker but I can while my phone is on loud speaker? I was on call with my friend and have put the phone on loud speaker and I opened youtube,The sound was audible to me but not to my friend.
How is it even possible when he can hear my voice but not the sound coming from the speakers

First of all, welcome to DSP SE.
In general most (if not all) smartphones nowadays have some kind of noise cancellation technology integrated. The simplest of all is to have two microphones, one which is supposed to record the noise and the other one record the voice with the noise (please keep in mind that this is an oversimplification of the problem). Then by subtracting the one from the other you end up with the voice without the noise.
Nowadays, there are some more sophisticated techniques to perform denoising such as echo-cancellation, active-noise-control schemes and many more, or even neural-network-based denoisers. In addition to that, the software running in the background may also perform some kind of speech recognition in order to get the speech components of the signal. Alternatively, some source-separation techniques could be used, but to be honest I am not sure how effective they could be.
All in all, the algorithm running in the background tries to separate your voice from the background (effectively, any other sound) and transmit only your voice. Most of the aforementioned techniques are adaptive filtering techniques, which means that they present a time-variable response. On top of that, they most probably will affect your voice. This is quite easy to experience yourself. Just try to move your phone around while you speak (or the person who you talk to should do that instead for you to hear the difference) or start and stop some sounds in the background abruptly.
I hope that this answer provides some insight into the general scheme of the process running in the background while you perform either phone or video calls. I know that this may not be a direct answer but it is a bit difficult to get an exact answer without providing the exact model and conditions (and to be honest I doubt you would still not be able to know exactly how it is done since this most probably constitutes the "magic" the company integrates into their products).


What are the different ways to find the envelope of a speech signal? I am new to signal processing.

I want to find the envelope of a given speech signal.

When I searched about it, I found about Mean Hilbert envelope function. I want to know what other functions/methods can do the same and why Mean Hilbert envelope is considered the best?

The notion of envelope is an ill-defined concept that describes the smooth upper and lower boundaries of more or less oscillating signals.

The most classical technique rectifies the signal (take the absolute value) and low-pass filters the result. Of course, results are not unique.

Other techniques identify local extrema and fit them with low-order functions like polynomials or splines. This is used for instance in some implementations of the Empirical Mode Decomposition.

The Hilbert transform sometimes offers a more solid alternative. Roughly speaking, suppose the signal can be written as the product of an envelope times an oscillation: $s(t) = e(t).o(t)$. Then, with other technical conditions, if the spectra of $e$ and $o$ are well separated, $H(s(t)) =e(t).H(o(t))$, where $H$ is the Hilbert transform. This is a loose transcription of the Bedrosian identity.

If $o(t)$ is a fine pass-band oscillations (the purest being a sine or a cosine), and $e(t)$ is lower in frequency, the Hilbert transform can separate them. This is why for speech, the signal is often decomposed first through a multirate filter bank.

All these techniques are often post-filtered in practice.

More recent approaches try to better define the regularity or smoothness of the envelope, using optimization techniques.


Transfer function of an Exponential system in Z domain Hi,
I am really confused with the system in the diagram.

The input-output relation of the system is given by $y[n]=\exp(x^2[n])$

I need to find the transfer function of this system $Y(z)$ in $z$-domain. My question is: is it really the simple exponential system I am missing out or something complex. If otherwise what would be the system transfer function.

This system doesn't have a transfer function, so there's no way to find it. Only linear time-invariant (LTI) systems can be described by a transfer function. In discrete-time, the $\mathcal{Z}$-transforms of input and output sequences of an LTI system are related by

$$Y(z)=H(z)X(z)\tag{1}$$

where $H(z)$ is the system's transfer function, which is the $\mathcal{Z}$-transform of its impulse response. Equation $(1)$ is another way to state that the output sequence is given by the convolution of the input sequence with the system's impulse response. If the system is not LTI there is no function $H(z)$ which relates $X(z)$ and $Y(z)$ as in Eq. $(1)$.

Your system is definitely non-linear so whoever asked you to find its transfer function is either clueless or wanted to fool you.


Opencv color segmentation of traffic signs: HSV ranges I am working on a project to detect and recognize signs in images from the LISA Traffic Sign Dataset. Since the images are very noisy (in regard to background) I am trying to segment the signs based  on their color using HSV and cv2.inRange().

The problem I am running into is my algorithm is segmenting far too much of the image. For example, when I try to segment the image based on different variations of yellow/orange for a "nice" image I get the following, without any false positives:

However, on images that are at a later point in the day, I receive the following:

I have played with the values for hours now, and cannot seem to find a happy medium that doesn't select both yellow and green in images. I was hoping the community could help provide some insight, since I've been scouring the web looking for examples of what I'm trying to accomplish for several days. I understand there are many examples of this being done in academic and scholarly paper/journals however I have no clue how to convert the complicated formulas they use into code.

Here is the relevant python code:

# the below is pairs of yellow and red respectively
hsv_color_pairs = (
    (np.array([21, 100, 75]), np.array([25, 255, 255])),
    (np.array([1, 75, 75]), np.array([9, 225, 225]))
)
...
# this code is in a for loop and loops over the above HSV color ranges
mask = cv2.inRange(hsv, colors[0], colors[1])
out = cv2.bitwise_and(im, im, mask=mask)
blur = cv2.blur(out, (5, 5), 0)
imgray = cv2.cvtColor(blur, cv2.COLOR_BGR2GRAY)
ret, thresh = cv2.threshold(imgray, 0, 255, cv2.THRESH_OTSU)

heir, contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

so to sum things up, my question is: what are some tips/tricks I can use to segment the images for the specific colors of the orange/yellow warning signs with regard to shadows/exposure without picking up other colors such as green.

Thanks all!

To those looking at this in the future, I realized I was being silly and forgot how HSV worked! To pick up the yellow more so than the green I needed to heighten the threshold of the saturation value within the range!

As a side note, I also was able to "normalize" the colors of the image a bit more by converting to Lab color space and applying CLAHE (for localized histogram equalization) and as well on the saturation values within the HSV color space. This allowed me to get the signs in ideal lighting conditions almost 100% of the time (roughly speaking of course) with very little false positives, but didn't really alleviate my problems on images with a redish/orangeish hue as seen in the second example image. But, as was mentioned by harshkn, this approach is very simple so I suppose I can only ask for so much success! Thanks to all who replied with suggestions!


Identify if a Specific Edge Pixel Belongs to Horizontal, Vertical or Diagonal Edge Line I have a code below that identifies the edge of the text image. Input image is:

Output display is:

I want to select random points near the edge border of the text. I was thinking to identify if the edge pixel is part of the horizontal or vertical line. If it  is part of the horizontal line, I can choose a vertical random point, else if the edge pixel is part of the vertical line, then I can choose a random horizontal point. How to it in Python? Aside from my logic, is there any better way to select the random points? Below is my initial code:

from PIL import Image
from matplotlib import pyplot as plt
import numpy as np
import cv2

image_counter = 1
path = 'sample/' + str(image_counter) + '.png'
image = Image.open(path).convert('L')
image_width, image_height = image.size
image = np.array(image)

sigma = .33
# compute the median of the single channel pixel intensities
v = np.median(image)
# apply automatic Canny edge detection using the computed median
lower = int(max(0, (1.0 - sigma) * v))
upper = int(min(255, (1.0 + sigma) * v))

#Edge Detection
edge_image = cv2.Canny(image,lower,upper)
plt.subplot(121),plt.imshow(image,cmap = 'gray')
plt.title('Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(edge_image,cmap = 'gray')
plt.title('Edge Image'), plt.xticks([]), plt.yticks([])
plt.show()

#Select random points outside text borders
#Select random points outside text borders

for x in range(image_width):
    for y in range(image_height):
        if (edge_image[x][y] == 255):
            #Check if the edge pixel is part of the horizontal or vertical line.
            #If part of the horizontal line, pick a random point vertically.
            #If part of the vertical line, pick a random point horizontally.
            # How to do it?

Desired output:

Not sure how you do it in python, but the idea can be as simple as the following:

After you select a random point on the edge pixel, calculate the horizontal gradient and the vertical gradient with this pixel at the centre of the gradient filter.

The width of the  kernel for the gradient filters will depend on the maximum width of the edges. For example if the edges are 3 pixel wide have a 5x5 gradient filter.

If the magnitude of vertical gradient is more than horizontal gradient then the pixel belongs to horizontal edge and vice versa for vertical edge.

Also for diagonal edge similarly you can have a kernel which calculated diagonal gradients, if diagonal gradients are not exactly 45 degrees, you might need to to increase the size of the kernel to get more angular resolution of gradients.

Hope that helps.


Condition for optimal value to be the average In many single parameter estimation process average is the optimal. Can anyone tell under what condition average is optimal? For example is average optimal if the data is iid, Gaussian or something else?

You don't need to assume anything about the data in order for the average to be optimal (in a least squares sense) for certain estimation problems. As an example, assume you want to estimate a constant value $x$ and you are given $N$ noisy measurements:

$$y_n=x+e_n,\quad n=0,1,\ldots,N-1\tag{1}$$

where $e_n$ is the error of the $n^{th}$ measurement. If our goal is to minimize the mean squared error

$$\frac{1}{N}\sum_{n=0}^{N-1}e_n^2=\frac{1}{N}\sum_{n=0}^{N-1}(y_n-x)^2\tag{2}$$

then the optimal estimate $\hat{x}$ is obtained by minimizing (2). Taking the derivative of the right-hand side of (2) with respect to $x$ and setting it to zero gives

$$\frac{2}{N}\sum_{n=0}^{N-1}(y_n-\hat{x})=0\tag{3}$$

From (3) the optimum estimate is obtained as

$$\hat{x}=\frac{1}{N}\sum_{n=0}^{N-1}y_n\tag{4}$$

which is simply the average of all observations.


Simulating a multipath channel I have a simple question, that might not be up to the level in this forum:

  In theory, a signal $s(t)$ that is received through a multipath channel reads:
  $$r(t) = \sum\limits_{k=1}^p \alpha_ks(t - \tau_k)$$
  where $r(t)$ is the received signal,  $\alpha_k$ is the gain of the $k^{th}$ path and $\tau_k$ is the delay of the $k^{th}$ path.

But

  In MATLAB, we have $s[n]$, sampled at some sampling period, $T_s$. To simulate $r[n]$, the only problem faced is that $\tau_k$ will not be a multiple of $T_s$. My question is "Is it a good idea to upsample the signal $s[n]$ up to a desired sampling period, say $NT_s$". If so, why should filtering be done after the upsampling process ?

The upsampling solution adds additional difficulties, as mentioned in some posts here. Please consider this solution..

The Discrete Fourier Transform (DFT) properties might support you. Take a look at the time shifting property. It states that $$\mathcal{F}\{x[k-\tau]\} = e^{-j\omega\tau}\mathcal{F}\{x\}$$ I implemented a possible use of the property in your problem in the following octave script. A picture with the result is also attached. A cos signal with 3[Hz] is sampled at 10[Hz], and shifted 33[ms] with the property. I hope it helps.

% sampled signal:
fs          = 10;
Ts          = 1/fs;
T           = 0:Ts:1;
nT          = length(T);
sig_sampled = 0.2 * cos(2*pi*3*T);

% just to plot, signal with high sample rate:
Tb          = 0:(1/1000):1;
sig_better  = 0.2 * cos(2*pi*3*Tb);

% delayed signal where tau is not a multiple of Ts:
% (using Fourier properties):
tau              = 0.033;
w                = 2*pi*(((-nT/2):((nT/2)-1))/nT)*fs;
fft_sig          = fft(sig_sampled);
fft_shift        = exp(-sqrt(-1)*tau*ifftshift(w)) .* fft_sig;
sig_delayed      = ifft(fft_shift);
sig_better_shift = 0.2 * cos(2*pi*3*(Tb-tau));

figure(1); clf;
plot(
    Tb, sig_better, 'k-',
    T,sig_sampled, 'k-x',
    Tb, sig_better_shift, 'r-',
    T, sig_delayed, 'r-x'
);
set(gca, 'ylim', 0.3*[-1,1])
grid on
legend('original signal', 'sampled', 'shifted signal', 'shifted sampled signal')
print('-djpg', '-r 200', 'fig.jpg')


How to invert FFT these magnitude equations to get FIR filters? I need filter designs that will reproduce or approximate each of the following:

1) ${\rm magnitude} = e^{-cf^{2}}$

2) ${\rm magnitude} = e^{-\frac{a+b\sqrt{gf}}{df}}$

Where: $a$, $b$, $c$, $d$ and $g$ are constants $>0$, and $f$ is frequency in Hz.

Or if it's easier or an equivalent problem, the multiple of the two filters as one combined filter would even be better:

3) ${\rm magnitude} = e^{-cf^{2} -\frac{a+b\sqrt{gf}}{df}}$

I have never done an inverse FFT or designed a FIR filter though so I'm looking for any basic guidance on how this might be done.

This is for physical modeling so the closer I can get the better. Accuracy is only important with audio frequencies (20-20,000 Hz).

I really need something that will work one way or another. Any guidance is appreciated.

As @RBJ commented the first one is a Gaussian filter. The Fourier Transform of a Gaussian is a Gaussian so this one can be easily created by sampling the impulse response given by the inverse Fourier Transform.

For case (1), an FIR filter with a magnitude response given by:

$$|H(f)| = e^{-cf^2}\tag{1}\label{1}$$

Given the Fourier Transform of a Gaussian:

$$ e^{-\pi t^2} \leftrightarrow e^{-\pi f^2}\tag{2}\label{2}$$

And using the Fourier Transform scaling property:

$$\frac{1}{|a|}x(t/a) \leftrightarrow X(af) \tag{3}\label{3}$$

We can rewrite $\ref{1}$ as:

$$X(af) = |H(af)| = e^{-\pi (af)^2}\tag{4}\label{4}$$

with $a = \sqrt{c/\pi}$

If we restrict $c$ to be a positive real number, then from $\ref{3}$ we get:

$$\mathscr{F}^{-1}\{ e^{\pi(af)^2}\}=\frac{1}{|a|}e^{-\pi(t/a)^2} = \sqrt{\frac{\pi}{c}}\left( e^{\frac{(\pi t)^2}{c}}\right)\tag{5}\label{5}$$

This is demonstrated with an increasing number of coefficients showing the excellent match that is achieved.

The MATLAB/Octave code used to create the filter coefficients was:

fs = 1;   # sampling rate
c = 100;
N = 27;   # number of taps, must be odd
t = [-(N-1)/2: (N-1)/2]/fs;   # time axis
coeff = sqrt(pi/c)*exp(-pi^2*t.^2/c); # sampled imp resp

Case (2) appears to be a high-pass filter which could be solved using the approach above by sampling the inverse Fourier Transform for $e^{-a/f}$ and $e^{-a/\sqrt{f}}$.

I couldn't easily proceed from this point and posted on the mathematics site for a possible solution:

https://math.stackexchange.com/questions/3666703/inverse-fourier-transform-for-e-a-f-and-e-a-sqrtf

A very simple but sub-optimum approach is to derive the filter coefficients from the inverse DFT. This is the frequency sampling approach to FIR design and provides an exact match at the frequency sampling locations, but results in significantly more ripple than optimized algorithms such as least-square and Parks-McLellan. The ripple can be significantly reduced by windowing the resulting coefficients, but for this application any windowing will significantly reduce the null at DC. However, in the spirit of "needing something that will work one way or another" the Frequency Sampling approach can be considered if a larger number of taps is feasible. Regardless of approach, the number of taps for an FIR filter is proportional to the inverse of the width of the transition band of the filter, so if the parameters are such that the filter is a tight DC-notch filter, a large number of taps with any traditional FIR approach will be necessary (an IIR filter may make a lot more sense, or a modification of Richard Lyons clever implementation of a linear=phase DC notch filter with CIC structure: https://www.dsprelated.com/showarticle/58.php ).

The example plots below show the result for a 275 and 501 tap FIR filter using the Frequency Sampling approach with the following example parameters as given in the code below to implement an FIR for Case (2):

fs = 1;
a = 500; b = 50; d = 50000; g = 100;
N = 275;
f = linspace(-fs/2, fs/2, N);
magC2 = e.^(-(a+b.*(sqrt(g*abs(f))))./(d*abs(f)));
coeffC2 = fftshift(ifft(ifftshift(magC2)));

The above plot appears to show an excellent match in both cases, however zooming in on the stop band shows that the frequency sampling approach only maintains a match (with the tighter parameters used) down to around -25 dB with 275 taps, while it maintains a match down to around -65 dB with 501 taps. Pursuing either an approach of sampling the computed impulse response from the inverse Fourier Transform, or a least squares algorithm to match the magnitude response will allow for a solution with significantly fewer taps. Also if the notch needed is actually not as tight as I have demonstrated here, the number of taps required will be significantly reduced as well.

Also observe how we can estimate the number of taps needed for any approach from the ideal in the plot above: if we want/need to maintain high accuracy out to -80 dB rejection, we observe for this case with the particular parameters chosen that the steepest transition the slope is approximately 65 dB over a frequency of 0.005 samples/cycle. From fred harris' rule of thumb this would be:

$$N \approx \frac{A}{22}\frac{1}{\Delta f} = (65/22)(1/0.005) = 590 $$

We can also see how increasingly challenging it will be to maintain a tight match at deeper rejections given the slope is continuously increasing.


Preamble Detection and Frequency Offset Estimation in radio packets I have generated a radio packet of data with a preamble header of 0x55555555 and crc32 for the data only.
The entire packet was separated in even and odd streams for I and Q. The I and Q were oversampled by a factor of 10.
The I and the Q were modulated into BPSK then then the I & Q were combined and mixed with the carrier and then transmitted. The transmission is random so I do not know when a packet starts or ends as the payload length is variable. The channel has multipath and noisy.

At the receiver I receive symbols with noise between the packets. I know there is frequency drift during transmission. I tried to frequency sync the packets, but I can not tell where the preamble starts from the collected symbols
If I demodulate the received symbols and search for the preamble I can not find it, because of the frequency drift.

What is a good approach to unlock this circular problem ?

Thanks.

For maximum sensitivity preamble detection with time referencing, consider using a barker code, or for more choices with longer lengths a PRN (pseudo-random noise) sequence, or for even more choices when more distinction is needed (CDMA) Gold codes are a possibility.  (Also not further detailed but have other advantages are Kasami codes and Zadoff-Chu codes).

All these codes offer the characteristic of having a strong autocorrelation at time offset = 0 and importantly a much smaller autocorrelation for any other time offsets. Such as that demonstrated in the graphic below for a PRN sequence using an LFSR (Linear Feedback Shift Register):

The codes provide some noise immunity through processing gain; under condition of white noise, the processing gain which is the improvement of SNR is given as $10Log_{10}(N)$ where N is the length of the sequence. This comes out from addition of correlated signal components results in the standard deviation of the signal component growing at rate N while the noise (if white, meaning uncorrelated from sample to sample) will have a standard deviation that grows at $\sqrt{N}$.  We detect the presence of the sequences using correlation, which is a multiply and sum:

$$Corr(\tau) = \Sigma r(t-\tau)x(t)$$

With the formula above representing a cross-correlation function of the received signal r(t) with the seequence of interest x(t), at various time offsets $\tau$.

There are several ways this is done in practice and depends on the detailed implementation approach but for shorter sequences this can be done with an FIR filter, with the coefficients of the filter as the code sequence upsampled to match the sample rate of the received signal (best to do this at the lowest rate possible of course).  For longer sequences it is often more practical to create a code generator with a programmable offset (to shift the start position in time) and simply multiply the received sequence by the code and accumulate the output. A threshold detection can indicate the successful presence of the code.

Below shows the result of a sliding correlation on a GPS signal, where the code repeats every 1 ms. (I would never use an FIR filter to correlate a GPS sequence which is 1023 chips long, but this shows what a sliding correlation would look like).

Other approaches involve doing the correlation in the frequency domain using this relationship:

$$Corr_{xy} = IFFT(FFT(x)FFT^*(y))$$

Which states that taking the inverse Fourier Transform of the complex conjugate product of the Fourier Transforms of x and y results in the circular correlation of x an y.

Further, a very important take-away is the sensitivity of any of these correlations to frequency offset, specifically that the magnitude goes down as a Sinc function with the first null at 1/T.

Please see this post directly where I detail that: GPS CA Signal Acquisition

This informs the frequency capture range based on the length of code used. If you are well inside this main lobe in frequency (I typically use 1/2 the null-null frequency range, but if the SNR is poor you may need to be more constrained). Typical acquisition approaches will either slowly ramp the frequency offset or step the frequency offset during correlation (The frequency is offset by using an IQ rotator on the received signal prior to correlation). Stepping the offset involves setting the frequency at an initial guess, correlating and if the acquisition threshold is not exceeded, stepping the frequency by 1/T (if using 1/2 the null-null frequency range) and repeating. If ramping, the ramp rate must be much less than 1/T without inducing an additional frequency offset loss.

The frequency offset is accurately estimated by comparing two successive correlations, since frequency is the change in phase versus the change in time ($d\phi/dt)$ by simply comparing the phase and time between two successive correlations, we can determine the frequency offset assuming we are in the frequency range of acquisition as detailed above (importantly these correlations are complex correlations with two real outputs I and Q such that carrier phase is determined).

Note that the frequency computation detailed at this link Obtaining phase and magnitude of 2 spectral components after ADC
using the computation

$$K = I[n]Q[n-1] - I[n-1]Q[n]$$  where K is propotional to frequency offset equally applies to the successive I and Q correlation outputs, referred to as I[n], Q[n] and I[n-1]Q[n-1].

Much more practical with modern processing and shorter preamble sequences would be a one shot frequency offset and delay estimation that I have detailed here:

GPS signal acquisition

That said, a robust preamble solution could involve sending multiple copies of one code to aid initial carrier and clock acquisition (and channel estimation!) followed by a final negative of the code sequence to mark the beginning of the data. (Just invert the code and the correlation will be shifted 180°). The trade space is over-head of the preamble versus acquisition sensitivity. It is typical to require more SNR during acquisition, especially if a fast acquisition is needed, hence the significance of correlation and processing gain.


How is the Tx/Rx Gain calculated in USRP? I know that USRPs are uncalibrated devices. In many exampls you can set the so called Tx or Rx gain (as in the picture below).

Now, I know that this does not correspond to the real power of the transmitter, nor is it a "gain" of the antenna - the most appropriate term here would be "transmitter channel gain", right?
So I would like to know how this transmitter channel gain process is implemented, how it is calculated, and how it works in practice (since a particular gain value doesn't correspond to a particular power level).

The term "gain" is pretty appropriate. It's a number describing the output power relative to some arbitrary base power.

So I would like to know how this transmitter channel gain process is implemented,

Depends on the USRP (or if you're using a modular model, on the daughterboard), but speaking generally: A mixture of selectable/bypassable amplifiers, variable gain amplifiers and adjustable attenuators.

how it is calculated,

Not calculated. Set by configuring the things mentioned above.

how it works in practice (since a particular gain value doesn't correspond to a particular power level).

It sets how much stronger or weaker your amplification is. Simple! You add 3 dB of gain, you get twice as much output power. Of course, that's approximate, and only within the linear range of all the signal chain.


Fastest Available Algorithm to Blur an Image (Low Pass Filter) Iam working with a camera that produces ugly artifacts:

by using ANY blur filter on the camera's output the visual quality improves drastically:

The above image was created using OpenCV's cv::medianBlur with a kernel size of 3.
I identified cv::medianBlur to be the fastest smooth/blur method in OpenCV.
However for my needs it is still too slow since it uses up to 80% of the whole processing time including encoding (ffmpeg MPEG4). I already tryed to use cv::UMAT but uploading each Image to the GPU and downloading the result again is taking even more time. So using OpenCL / Cuda isn't an option either!?

Therefore Iam looking for the fastest way to blur/smooth an image. However there are so many algorithms in so many librarys(OpenCV/IPP/swscale) to implement and test them all would take to much time. So do you have any suggestions which algorithms I can take a look at, or could offer a really good performance?

Here are some test results for 75.000 iterations of a 640x360 image:

+------------------------------+--------+----------+
|          Algorithm           | Kernel | Time(ms) |
+------------------------------+--------+----------+
| cv::medianBlur               | 3x3    | 18492    |
| cv::medianBlur ocl           | 3x3    | 54596    |
| ippiFilterMedianCross_8u_C3R | 3x3    | 15755    |
| cv::blur                     | 3x3    | >100000  |
| cv::GaussianBlur             | 3x3    | >100000  |
| cv::filter2d                 | 3x3    | >100000  |
+------------------------------+--------+----------+

The fastest blur would be Box Blur.
You can implement it using Running Sum.
I think Intel FilterBoxBorder works in that manner.

If you'd like you can do a few passes of it to approximate the Gaussian Blur.

You can also use IIR Filter Coefficients to blur the image quite easily.

You may have a look at my project Fast Gaussian Blur.


Modulate Bipolar(?) signal I want to modulate this signal if possible with an Arduino, but have no idea what circuit can help me here. Let's ignore the timing and the measurements from the image for a moment, my problem is more about the alternating shape of the signal.
I'm a beginner here and learning a lot at the moment. I don't wan't any complete circuits or code, but would be grateful if u could push me in the right direction. Thanks

For a signal that rides at VSS/2, goes down to 0V and up to VSS, just use a pair of GPIO pins, each with a $10\mathrm{k\Omega}$ resistor, with the far ends connected together.  When both pins are high you'll get VSS, when both are low you'll get 0V, when one is high and the other low you'll get VSS/2.
Then sequence the pins to get the pattern you want.  If you need something that's faster and has less jitter than you can accomplish by bit-banging the pins, use a pair of timer outputs and sequence them in hardware (which you may not be able to do with the ATMega found in most Arduinos, and which the Arduino hand-holding software may not allow you to do easily -- it's your job to figure that part out).
To get your 2V, 3V, 4V pattern, you'll need up to four resistors, connected together in the middle, and to VSS, GPIO_n, GPIO_{n+1}, and ground on the other ends.  Finding their value takes some algebra, but it's straightforward if you know basic circuits.


Why is LTE either FDD or TDD? In LTE systems, up- and downlink are duplexed either using FDD or TDD. Why has this not been designed to utilize the same multiplexing technique as used for the downlink (OFDMA)? Resource blocks could be dynamically assigned for up- and downlink then, making everything much mor flexible.
I know that for uplink sc-fdma is utilized, to minimize PAPR, is this maybe the only reason?

You are mixing two things modulation and multiple access. SC-FDMA has been chosen because it is battery-friendly modulation. SC-FDMA is used in both FDD and TDD modes.

There are several reasons why both FDD and TDD are defined.

Historically, people want to utilize efficiently the precious spectrum. FDD needs two bands while TDD needs only one.

Economically, LTE is a broadband technology, i.e. it is designed for high speed human data connection, and high speed human data connection favors downlink. For example in China, the dominant (and at the beginning of commercial LTE, the only) market of TDD LTE, Youtube video streaming takes a large part. In this case, reserving a band for UL is kind of wasteful. Moreover, LTE was commercialized in the middle of the Great Recession and operators were very skeptical in CAPEX spending.

Protocally, FDD and TDD are very different in term of multiplexing devices (UE), especially in HARQ-ACK reporting, UE scheduling etc. If you look at the procedure standards of LTE, just PHY layer for example, each mode requires very complicated scheduling. In mixing two modes in one cell, the level of complexity is huge for eNB, with no relevant benefit. However, mixing two modes between two or more cells are considerable and has been standardized since Release 11 in the Carrier Aggregation scope.

Technically, TDD supports better channel information (CSI) reporting thanks to the reciprocity of channel. Thus it is favored for beamforming, which requires accurate channel estimation, and a very good candidate for massive MIMO.

Finally, yes you can mix two modes if you want, but the efforts would not worth the benefit you get.


Creating a PAM function with Matlab So I need to make a function that will be called as function [t, pulse] = PAM (array). My function should have a result like this when I use array = [1 0.3 0.8 0.6], which array is the width of my digital pulsing signal.

Also please not a complete answer, I think it's better to think alone, but I'm in a point that I don't know what to do.

Solved it!
function [t, pulse] = PAM (array)
  syms t
  u(t) = heaviside(t);
  pulse = array(1)*u(t);

  for i = [2:length(array)]
    if array(i) > array(i - 1)
      pulse += abs(array(i-1)-array(i)).*u(t - (i-1));
    else
      pulse -= (array(i-1)-array(i)).*u(t-(i-1));
    endif
  endfor
endfunction


Exponential decaying Swept / chirp signal I'm trying to create a chirp/swept signal that has it's frequency going down exponentially using matlab/octave.  I can't use the built in chirp function because I'm trying to create a chirp/swept signal based on an equation/formula.  I tried creating some code below but as you see from the spectrogram the frequency does not end at the correct frequency, and the frequency doesn't go down exponentially.

clear all,clc
freq1=20; %start freq
freq2=200; %end freq
fs=44100;
dur=1; %duration of signal in seconds

t = linspace(0,2*pi,fs*dur);

f=freq1:(freq2-freq1)/length(t):freq2-(freq2-freq1)/length(t);
%20:(200-20)/lenght(t) :200-(200-20)/length(t)

data=sin(f.*t); %build signal
data=(data/max(abs(data))*.8); %normalize signal
wavwrite([data'] ,fs,32,strcat('/tmp/del.wav')); %export file
plot(t,data)

Here's a spectrogram plot

My goal:

1) create a signal based on a given equation/formula where the frequency is exponentially decaying (based on the formula 1/exp)

2) be able to type in the starting and ending frequency along with it's duration in seconds example
startfreq=20; endfreq=200; durseconds=2;

The spectrogram plot should look something like this if the start frequency is 200hz ending at 20hz.

Any ideas?

Frequency is the derivative of phase, so you need to start with the expression for frequency, which in your case is
$f=200\exp(-\alpha*t)$

At time $t=0$ this will have the frequency of 200, but we need to solve for $\alpha$ so that the frequency at time $t=2$ is 20.  Solving for $\alpha$ gives $\alpha=1.1513$.
To get the correct expression for the phase you need to integrate the above expression for the phase i.e. solve

$f =\frac{1}{2\pi}\frac{d\phi}{dt}$

The following matlab code generates this signal.

fs=500; %sampling frequency
t=(0:2*fs)/fs; %Time vector
alpha=1.1513;
sig = exp(-j*2*pi*200/alpha*exp(-alpha*t));

Note that just because the instantaneous frequency behaves this way doesn't mean the spectrogram will follow it exactly. It will get better if you have a higher time-bandwidth product.


Extreme quantization of audio I would like to take an audio sample and reduce it to a very limited set of numbers, with one value for amplitude and another for length.

Example:

I want to find the dimensions of the red boxes, which represent areas of similar sound intensity. Ideally I'd like this as a python list in a format similar to:

[[2,10],[5,3],[2,8], ...]

(where [x,y] = [length,intensity])

During creation I'd like to specify resolution in some way, so I can take one sample and have a list of 100 items or 10 items or whatever.

Being able to specify y resolution would also be useful, so \$intensity can be 1-10 or 1-100 or so on. As I write that I realise that I can probably only specify one of those, if \$intensity is bounded as an integer 1-10 then the number of list items depends on that, I think?

I'm pretty familiar with using python for various things, audio/signal processing is totally out of my area of expertise so I'm a bit lost and any pointers on where to start would be much appreciated.

The context is I'd like to use that list as an input in a FreeCAD or Blender script which will generate various shapes and designs from audio samples, ultimately to be 3D printed.

So, you describe what you want to do relatively explicitly:

You care about the amplitude, not the sign of the signal.
you need to detect periods where that doesn't change much.
You want to quantize the result.
You want to Run-Length Encode (RLE) the result

So, to achieve each point:

abssignal = numpy.abs(signal)
Apply a low-pass filter. Simply design one using scipy.signal.
Quantization as ${0,\ldots,N}$: Divide by max(filtered_signal), multiplying by $N$ and then casting to integer (numpy.astype).
Apply RLE


Energy using Hanning window I am a newbie in signal processing. I am testing the Hanning window on a signal. Obviously doing this window also means cut some energy of the signal, and when I do the DFT, the 0 frequency is not equivalent anymore to the average of the time history. Do you usually solve in some way this problem?
Can someone suggest papers or books where noise-reducing strategies for the FT of a signal are investigated? I just don't know how to start, an online search didn't give me more than basic definition (so please don't send me the Wikipedia definition of the tapering window)

Best

I suspect you are talking about the Hann window, after Julius von Hann, and Blackman & Tukey. Of course, a non-constant window may induce energy variations, either used as a weighted average, or a convolutional filter. This can be compensated in several ways:

by normalizing the window amplitudes
by separating tipme-portions, or signal components, like removing the mean or low-frequencies before doing the processing
by modelling the signal and applying model-based corrections.

More would be needed on your signal to provide more precise answers. R. Lyons on a diagram on page xviii of his book Understanding Digital Signal Processing highlights the importance of windows.


Why resampling change this FFT output? I have been studying the following code but still unable to explain its output. I really hope someone could shed some light.

Let me begin with the baseline. Suppose I have a 20 seconds sinusoid with frequency components of [1,2,3,4....500Hz]. The signal has Fs=1024Hz, so it's FFT output looks as follows. The FFT amplitude is also scaled to preserve Parseval's Theorem. The frequency resolution is Fs/nfft = 1.0 Hz per bin in this case, as the result the FFT output smeared itself and make the output looks like a flat line, which is expected.

# Python code

import numpy as np
import matplotlib.pyplot as plt

dur = 20
fs = 1024
t1 = np.arange(0, dur, 1/fs)
x1 = np.zeros_like(t1)
for f in np.arange(1, 500):
    x1 += np.sin(2*np.pi*f*t1)

plt.figure(figsize=(8,4))
plt.title('Signal')
plt.plot(t1[:100], x1[:100])
plt.xlabel('sec')
plt.grid()
plt.show()

nfft = 1024
df = fs / nfft
X1 = np.fft.fft(x1, nfft)[:nfft//2]

# Scale FFT amplitude to preserve Parseval's Theorem
X1 = X1 * np.sqrt(2 / nfft)
f1 = np.arange(0, nfft//2) * df

plt.figure(figsize=(8,4))
plt.title(f'Baseline FFT, Freq Resolution={df}Hz')
plt.plot(f1, np.abs(X1))
plt.grid()
plt.xlim([0, nfft/2])
plt.xlabel('Hz')
plt.show()

p1 = np.sum(x1[:nfft]**2)
p2 = np.sum(np.abs(X1)**2)
print('Parseval check:', np.allclose(p1, p2))

Here comes the tricky part. The signal is resampled in the following way, x2 = x1[::10] (downsampling by a factor of 10?). If that's really downsampling then new effective Fs should become 1024/10 = 102.4Hz, then the Nyquist frequency becomes $102.4/2 = 51.2Hz$, so the FFT output should only contain frequency components up to 51Hz. At least that's from the theoretical perspective. However the following output suggests otherwise. I can see the frequency resolution has improved significantly from 1.0Hz to 0.1Hz, but then the output only contain frequency components up to 12Hz instead of 51Hz. Why is that happen? What do I miss? Can anyone shed some light? Thanks

step = 10
x2 = x1[::step]

plt.figure(figsize=(8,4))
plt.title('Downsampled Signal')
plt.plot(t1[:100], x2[:100])
plt.xlabel('sec')
plt.grid()
plt.show()

nfft = 1024
fs2 = fs / step
df2 = fs2 / nfft
X2 = np.fft.fft(x2, nfft)[:nfft//2]

# To preserve Parseval's Theorem
X2 = X2 * np.sqrt(2/nfft)
f2 = np.arange(0, nfft//2) * df2

plt.figure(figsize=(8,4))
plt.title(f'Downsampled FFT, Freq Resolution={df2}Hz')
plt.plot(f2, np.abs(X2))
plt.grid()
plt.xlim([0, 50])
plt.show()

p1 = np.sum(x2[:nfft]**2)
p2 = np.sum(np.abs(X2)**2)
print('Parseval check:', np.allclose(p1, p2))

I think there are a couple things going on here.

First, the flatness of your baseline plot is a little deceiving because the tones that you have added (1 - 500 Hz) are all almost exact multiples of the fundamental frequency of the fft (fs/nfft = 1.024 Hz/bin).

If you increased nfft by a factor of 10, for example, the finer frequency resolution of the zero-padded FFT would show something similar to your New FFT Plot.  Here it is with nfft = 10240:

Because the tones are centered well on the FFT bins, there is no spectral leakage.  If you were to change the tones so that they were for example at multiples of 0.75 Hz, you'd see something like this:

It looks like a mess, but what's really happening is that each tones is getting smeared across several bins and the tones are constructively and destructively interfering.  Here's a single tone at 0.75 Hz with nfft=10240:

This spectral leakage occurs because the FFT assumes that its input is repeated forever, so if you have a sinusoid whose period is a multiple of the fundamental frequency of the FFT everything looks as you'd expect, but if not then there will be a discontinuity between the first sample and the last sample.  This effect can be mitigated using window functions: https://en.wikipedia.org/wiki/Window_function.

However, the main issue that you're seeing is the result of aliasing.  If you want to avoid aliasing, you need to low-pass filter your signal before downsampling.  Otherwise, any frequency content in the original signal that is above the Nyquist rate of your down-sampled signal will alias.  In your case, the sampling rate after down-sampling is 1024/10 = 102.4 Hz, so all of the tones you have above 51.2 Hz in your original signal will alias to frequencies below 51.2 Hz and interfere with the true signals at those frequencies.

Changing your down-sampling factor to 8 makes the effect easier to demonstrate because the aliased tones will land right on top of other true tones.  If your original signal only has tones below the 64 Hz, all tones are present in the down-sampled signal:

for f in np.arange(1, 64):
...
step = 8
x2 = x1[::step]   # down-sampling by a factor of 10??
fs2 = fs/step

Now let's add a tone at 65 Hz:

for f in np.arange(1, 66):

There is no longer any tone at 63 Hz or 65 Hz.  The 65 Hz tone has aliased and destructively interfered with the true tone at 63 Hz to cancel it out:

Now, let's change the phase of the original tone by 180 degrees:

for f in np.arange(1, 66):
x1 += np.sin(2*np.pi*f*t1 + np.pi/2)

There's still no tone at 65 Hz, but the one at 65 Hz is now constructively interfering, doubling the magnitude in the 63 Hz bin.

So if you extend this to 500 Hz, you can start to see what's going on.  The fix is to to apply a low-pass filter before you down-sample.  The combined steps of low-pass filtering and then down-sampling is called decimation.  Here's the scipy function to do it: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.decimate.html.


Channel Estimation / Equalization - Estimate Channel Inverse Using White Noise Statistics Only Given the system defined in the following figure:

We have a system $ G \left( f \right) $ which is unknown yet can be defined by $ {N}_{p} $ poles and $ {N}_{z} $ zeros.

The signal $ x \left( t \right) $ is an IID White Noise.
This signal samples can not be accessed.

The signal $ y \left( t \right) $ is the colored noise (Filtered).
This signal samples are available.

The filter $ H \left( f \right) $ is adaptive (Time changing).
Its model, which is to be determined, is given by $ {N}_{z} $ poles and $ {N}_{p} $ zeros.

The signal $ z \left( t \right) $ is the result of the $ H \left( f \right) $ Filter.
This signal samples are available.

How can one set the parameters (Poles and Zeros) of $ H \left( f \right) $ such that $ z \left( t \right) $ is White Noise?
This means how to adaptively set the parameters of $ H \left( f \right) $ to be the inverse filter of $ G \left( f \right) $.

Usually those kind of problems are solved using the Least Mean Square Filter yet since the samples of $ x \left( t \right) $ aren't accessible it is not the classic model.

Is there any other way to set the parameters of $ H \left( f \right) $ s.t. its output, $ z \left( t \right) $, is white noise without requiring access to $ x \left( t \right) $?

A Regular Random Process is is the the result of White Noise Going through a Minimum Phase LTI System.

A Non Perfectly Predictable Random Process can be defined as (See Wold Theorem / Wold Decomposition):

$$ x \left[ n \right] = \sum_{k = 0}^{\infty} {g}_{k} u \left[ n - k \right] + z \left[ n \right] $$

Where $ {\left\{ {g}_{k} \right\}}_{k = 1}^{\infty} $ is a Minimum Phase LTI System, $ u \left[ n \right] $ is White Noise and $ z \left[ n \right] $ is a Perfectly Predictable Process.

Since the questions deals with colored noise, $ z \left[ n \right] = 0 $.

Hence the process can be defined as following:

The optimal linear predictor is given by:

$$ \hat{x} \left[ n \mid n - 1 \right] = \sum_{k = 0}^{\infty} {h}_{k} x \left[ n - 1 - k \right] $$

Where the filter $ H \left( f \right) $ is given by:

$$ H \left( f \right) = 1 - \frac{1}{G \left( f \right)} $$

Hence, given the Linear Prediction Filter one could define the Whitening Filter as:

$$ W \left( f \right) = 1 - H \left( f \right) = \frac{1}{G \left( f \right)} $$

The Linear Predictor Filter can be approximated by the LMS Filter:

Which now, the Whitening Filter can be easily defined and approximated using the LMS Framework.


Solving a Linear Mean Square Estimation the Easy Way I have an exercise which is quite trivial. However I got stuck and I'm not sure if this the end-result. I assume there has to be a way to get this result much quicker.

Given are two randomly distributed variables y and n with mean 0 and variance 1. Further we know that $E\{yn\} = 0.5$ We measure $$x=y+n$$

We look for linear-mean-square estimation of y as a function of x.

From the given facts we now x has zero mean two. The variance has to be $$var(x) = var(y)+var(n)+2\cdot cov(yn)$$
because y and n are correlated. Therefor $var(x)$ becomes $$1+1+2\cdot 0.5 = 3$$

y becomes $y = n-x$ and $\hat{y}$ should have the form $\hat{y}= ax+b$.
For a and b I already have the general formula:
$$a=\frac{E\{xy\}-m_xm_y}{\sigma_x^2}$$
$$b=\frac{E\{x^2\}m_y-E\{xy\}m_x}{\sigma_x^2}.$$
In my case y should be $n-x$ if i got this right. But solving the problem with those formulas for a and b leaves me with terms like $E\{xn\}$. How can I know what this mean should be? I can't assume they are uncorrelated, can I? Probably I'm running completely in the wrong direction, and the solution is much more obvious.

EDIT 1: From Fat32's input I get for a:
$$a = \frac{E\{x(n-x)\}}{\sigma_x^2}=\frac{3-1.5}{3}=1/2$$
and b $$b=\frac{E\{x^2\}m_y-E\{x(n-x)\}m_x}{\sigma_x^2}=\frac{3-1.5}{3}=\frac{1}{2}.$$
The solution would therefor be $y=-\frac{1}{2}x+\frac{1}{2}=\frac{1}{2}(1-x)$ Not sure if this is true. Have to test it with random samples.

EDIT 2:
I did a test with matlab:

N = 10000;             %// Number of samples in each vector
M = randn(N, 2);

R = [1 0.5; 0.5 1]; %// correlation matrix
M = M * chol(R); %used to calculated depended random variables

n = M(:, 1);
y = M(:, 2);

x = y+n;

y_hat = -0.5*x+0.5;

mean(y-y_hat)

y_hat is not even close to the real y. Has not even the same mean. I don't get it. I'm making definitely some mistakes here.

EDIT 3:
Found another formula which uses the a and b. Inserted the linear leas squares solution becomes $$\hat{y}=\rho_{xy}\frac{\sigma_y}{\sigma_x}(x-m_x)+m_y.$$ When I insert my values I get:
$$\hat{y}=\rho_{xy}\frac{1}{\sqrt{3}}x.$$
Rho is $\frac{1.5}{\sqrt(3)}$ and so $\hat{y}$ becomes 0.5x as @Fat32 pointed out. The error above was that b is zero because $m_x$ and $m_y$ are zero.

So in your case doesn't the relation $x = n+y$ help ?

I mean, assuming your derivation for the mean square estimtor is right, then to compute $E\{xn\}$ you would look for $E\{ (y+n)n\}$ and using properties of $x$ and $n$ you would get

$$E\{xn\} = E\{(y+n)n\} = E\{yn\} + E\{n^2\} = 0.5 + 1 = 1.5 $$


Why is a linear phase important? If symmetry conditions are met, FIR filters have a linear phase. This is not true for IIR filters.

However, for what applications is it bad to apply filters that do not have this property and what would be the negative effect?

Let me add the following graphic to the great answers already given, with the intention of a specific and clear answer to the question posed. The other answers detail what linear phase is, this details why it is important in one graphic:

When a filter has linear phase, then all the frequencies within that signal will be delayed the same amount in time (as described mathematically in Fat32's answer).
Any signal can be decomposed (via Fourier Series) into separate frequency components. When the signal gets delayed through any channel (such as a filter), as long as all of those frequency components get delayed the same amount, the same signal (signal of interest, within the passband of the channel) will be recreated after the delay.
Consider a square wave, which through the Fourier Series Expansion is shown to be made up of an infinite number of odd harmonic frequencies.
In the graphic above I show the summation of the first three components. If these components are all delayed the same amount, the waveform of interest is intact when these components are summed. However, significant group delay distortion will result if each frequency component gets delayed a different amount in time.
The following may help give additional intuitive insight for those with some RF or analog background.
Consider an ideal lossless broadband delay line (such as approximated by a length of coaxial cable), which can pass wideband signals without distortion.
The transfer function of such a cable is shown in the graphic below, having a magnitude of 1 for all frequencies (given it is lossless) and a phase negatively increasing in direct linear proportion to frequency. The longer the cable, the steeper the slope of the phase, but in all cases "linear phase".
This makes sense; the phase delay of 1 Hz signal passing through a cable with a 1 second delay will be 360°, while a 2 Hz signal with the same delay will be 720°, etc...
Bringing this back to the digital world, $z^{-1}$ is the z-transform of a 1 sample delay (therefore a delay line), with a similar frequency response to what is shown, just in terms of H(z); a constant magnitude = 1 and a phase that goes linearly from $0$ to $-2\pi$ from f = 0 Hz to f = fs (the sampling rate).

The simplest mathematical explanation is that the a phase that is linear with frequency and a constant delay are Fourier Transform pairs. This is the shift property of the Fourier Transform. A constant time delay in time of $\tau$ seconds results in a linear phase in frequency $-\omega \tau$, where $\omega$ is the angular frequency axis in radians/sec:
$$\mathscr{F}\{g(t-\tau)\} = \int_{-\infty}^{\infty}g(t-\tau)e^{j\omega t}dt$$
$$u = t - \tau$$
$$\mathscr{F}\{g(u)\} = \int_{-\infty}^{\infty}g(u)e^{-j\omega (u+\tau)}du$$
$$ = e^{-j\omega \tau}\int_{-\infty}^{\infty}g(u)e^{-j\omega u}du$$
$$ = e^{-j\omega \tau}G(j\omega)$$


Allpass Filters - Causal and Stable So I have been learning about how to test systems for causality and stability but I am confused about the implications on their unit circle representation.

Would it be safe to say that a causal and stable allpass filter does not have a zero inside the unit circle?

You're right, a causal and stable discrete-time system must have all its poles inside the unit circle. Since an allpass system with a zero $z_0$ inside the unit circle ($|z_0|<1$) must have a pole at $1/z_0^*$, i.e., outside the unit circle, it cannot be causal and stable.


What is causing these spikes in the output of the low pass filter? I am low pass filtering the PCM of an arbitrary Bell 103 (300 bps) signal. However, the output of the low pass filter shows random spikes in amplitude.

I have tried to increase the bandwidth, but no luck.
Can someone explain why these spikes arise? My goal is to do some sort of envelope detection on this signal afterwards, to decode the 1's and 0's. For example by using a rectifier first.

The filter coefficients:

coeffs = [0.265842424474100 0.243276902464240 0.271096721997119 0.243276902464240 0.265842424474100]

The Matlab script used:

close all;
clear all;

% Generate a random bit sequence as input
bitsIn = randi([0 1], 16, 1)';

% Generate the amplitude of the FSK signal
amplitude = randi([7000 10000], 1, 1);

% Generate a randomized length of randomized startup noise
% x = randi([-20 20], 1, randi([0 3000], 1, 1));
x = []; % TODO

% FSK module the input bit stream
for i = 1:length(bitsIn)
    if (bitsIn(i) == 1)
        freq = 1270;
    else
        freq = 1070;
    end

    % Append the new PCM data for the new bit, with random phase
    phase = (2 * rand(1)) * pi;
    x = [x (amplitude * sin((2 * pi * freq * (0:26) / 8000) + phase))];
end

% Lowpass filter the resulting FSK samples
d = fdesign.lowpass('N,Fp,Fst,Ast',4,1070,1270,10,8000);
Hd = design(d,'equiripple');
y = filter(Hd,x);

figure();
plot(x);
title('Input data');
figure();
plot(y);
title('Low pass output');
```

I don't have MatLab with the DSP toolbox, but you're calling filter wrong.  Use:

y = filter(Hd, [1], x);

to set the feedback taps properly to $[1]$.

You should also inspect your filter design by using

freqz(Hd, [1], 512, 8000);

Although my hunch is that amplitude spikes in the output are directly caused by the filtering of the phase discontinuities you intentionally added when performing FM modulation.

Really, no modern equipment uses frequency switching to perform FSK.  FSK emitted from modern equipment is almost always continuous phase.  Are you really dealing with an ancient piece of Bell 103 modem hardware?

An example Octave script that generates the complex I/Q samples of continuous phase FSK is here:

https://dsp.stackexchange.com/a/53834/28112

if you need an example of how to generate CPFSK.


fft(x) different from fft(ifft(fft(x))? I normalized the spectrum of a time series and windowed the spectrum, then something strange happened. The code below can run in MATLAB or Octave.

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);
k = n/10;
w = sin(linspace(0,pi/2,k));
y1(1:k) = y1(1:k) .* w;
y1(m:-1:m-k+1) = y1(m:-1:m-k+1) .* w;
y1(n-m+3:n) = y1(m-1:-1:2);
x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');

As shown in subplot 3, the red line is the difference between x1 and x2, meaning that they are exactly the same, but as shown in subplot 4, their spectrum (black and red, respectively) is different.
So, It's weird FFT of IFFT of a spectrum is different from the spectrum itself.

The problem appears to be that your windowing function doesn't preserve complex conjugate symmetry. So ifft(y1) has a significant imaginary part. By discarding this through the real() operation you induce a significant error which results in the discrepancy. To verify try

z = ifft(y1); plot(imag(z));


Confused about Fourier transform output I am new to signal processing, so I did a fourier transform on a frame of 200 samples with NFFT = 512 and then I took the absolute value :

I am trying to understand these values on the output, are these frequencies ?
I tried to plot them, here is the result :

Does it mean that my frame contains 60khz frequencies, I am a little confused about this output.

The frequencies at FFT output are relative to the input sampling rate fs and the FFT size N. The $k^{th}$ index corresponds to a frequency $k*fs/ N$. Pay attention that values above $N/2$ correspond to negative frequencies


Do Nyquist samples have the same sum as the integral of the signal? Assume I sample a signal according to the Nyquist criterion. Then I perform a simple summation / integral over a linearly interpolated signal. Is this equivalent to the integral over the continuous signal, i.e.?

$$\sum_{n=-\infty}^{\infty}x(nT) T = \int_{-\infty}^{\infty}x(t)\text{d}t$$

I hope my representation is correct for what I mean. My question relates to whether a reconstruction, e.g. Whittaker–Shannon interpolation, is necessary   in practical application. I would also appreciate any literature pointers for discussions on how this might affect practical cases with finite time and imperfect filters.

The equation in your question is correct. If we assume that $x(t)$ is an ideally band-limited signal with bandwidth $B$, then, according to the sampling theorem, we have

$$x(t)=\sum_{n=-\infty}^{\infty}x(nT)\frac{\sin\left[\pi (t-nT)/T\right]}{\pi (t-nT)/T}\tag{1}$$

where the sampling rate $1/T$ satisfies $1/T>2B$.

If we assume that the integral over $x(t)$ exists, using $(1)$ we obtain

$$\begin{align}\int_{-\infty}^{\infty}x(t)dt&=\int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}x(nT)\frac{\sin\left[\pi (t-nT)/T\right]}{\pi (t-nT)/T}dt\\&=\sum_{n=-\infty}^{\infty}x(nT)\int_{-\infty}^{\infty}\frac{\sin\left[\pi (t-nT)/T\right]}{\pi (t-nT)/T}dt\\&=\int_{-\infty}^{\infty}\frac{\sin(\pi t/T)}{\pi t/T}dt\cdot \sum_{n=-\infty}^{\infty}x(nT)\\&=T\sum_{n=-\infty}^{\infty}x(nT)\tag{2}\end{align}$$

where the last equality follows from the fact that the integral in the one but last line just equals the DC value of an ideal low pass filter with gain $T$.


Power Spectral Density as a single number confusion I'm trying to recreate the results of a machine learning applied to the DSP classification in the article: link.
I have a signal (activity measurements from a smartwatch) per patient, so about 30 signals in total, quite long. For machine learning algorithm authors take traditional approach: extract statistical values about the signal, create columns (1 column per feature, 1 row per patient/signal) and plug this into a classifier. I'm confused about the spectral density used as feature.
Questions:

If I take PSD of a signal, I get another signal (time series) with values, not a single value, am I right? Therefore I would get not 1 column with spectral density feature per patient, but a full signal.

Suppose the article is imprecise (it is in other places) and authors made some "mental shortcut" while writing this. What else could they mean by "spectral density" feature? It should be 1 column, so 1 number extracted from each signal from patient. They mention it is in frequency domain, so PSD would be right, but it would not produce a single number.

Question 1: Yes you are correct, the power spectral density is the power distribution per unit frequency so is a continuous function of frequency.
Question 2: The single number as given can be an estimate of total power.
What they gave is completely incorrect starting with the formula as given:
$$P = lim_{T \rightarrow \infty}\frac{1}{T}\int_0^T|x(k)|^2dt$$
This is an attempt to provide a formula for the total power in the signal, as the time average of the energy of the signal (but where is t in the function for x?). I assume $k$ is a typo, but if it refers to an FFT bin, then as $T \rightarrow \infty$ discrete $x(k)$ would become a continuous function of time as $x(t)$ so that the above would properly represent the total power as:
$$P = lim_{T \rightarrow \infty}\frac{1}{T}\int_0^T|x(t)|^2dt$$
Note that the absolute value squared of the FFT bins, $|x(k)|^2$ can be used as an estimate of the power spectral density (see links below for conditions on this), with proper scaling to any units desired as the the relative power in each bin.
For further details on estimating PSD from the DFT, see:
Power spectral density vs. FFT bin magnitude
What is the difference between PSD and squared magnitude of frequency spectrum?


RF-Chain Signal Delay for Sensor Switching Lets have a RF-chain as above with system bandwith from HP corner to the green line say 1 MHz.
The signal accumulates delay as it passes through this analog chain, and due to it being non-liner filter chain, the group delay is likely as shown. Now at the end of LP1 and LP2 we have some form of detectors, which taken in the signal and do some type of functional detection on it. However the most important part of the detection is the critical control of switches s1 and s2 to determine when the signal is actually valid ( coming out of the two LFs ?).
So, the question would be what is a valid signal out of the LPF1 , since the input can not be a fixed single freq. component, otherwise we could just use the phase delay. Here the bandwidth is large upto 1 Mhz.
So essentially we have to determine the optimum read-out freq or some statistical values to find the optimum switching times.
Can you please suggest how can one actually go on  to do it ?
Without any new hardware design . Just lets try signal processing and statistical analysis
Thanks

The OP clarified in comments under the question that the intention is to compensate for group delay distortion introduced in the hardware. The typical approach to optimally compensate for this (in a least squared sense) with processing alone is to use the Wiener-Hopf equations to determine the coefficients of an equalizer that can be implemented with an FIR filter (meaning a difference equation with only feed-forward terms). I detail the full approach of doing this in these posts so will provide the links below, but to bottom line the process, the channel to be equalized (the receiver) is "sounded" with a known waveform that is spectrally rich (pseudo-random noise or frequency chirps are great choices as they also offer high average power or SNR; an impulse is a poor choice since it is a challenge to do that with high SNR), and with that and the received signal after the channel the reverse deconvolution can be computed in a least squared sense to determine the effective inverse channel, but importantly applicable to mixed phase systems (systems with both leading and trailing echoes) which on their own can't be inverted due to having zeros in the right-half plane. This approach would be ideally suited for distortions introduced in hardware implementations that are not changing with time (within an acceptable tolerance) since the channel can be sounded once and fixed coefficients determined relatively easily with pre-processing that can then be used without further modification, in contrast to the iterative algorithms that are needed to provide an adaptive solution for time-changing channels.
Details of the Wiener-Hopf Equations and shows application to determine the transfer function of the channel:
Compensating Loudspeaker frequency response in an audio signal
Equalizer Implementation Example by swapping Tx and Rx in the previous case, we can instead solve for the causal equalizer for a mixed phase system instead of the channel itself:
How determine the delay in my signal practically


Split and Merge image segmentation algorithm in OpenCV I asked google and I've tried to find something in the API of OpenCV, but hadn't any good hints on split and merge image segmentation method.
I want to test it, if it's working better than my proposed method which I posted a while ago here on SP.
Do you know if there exists a free library or OpenCV extension for this segmentation method?
Any good hints for writing it fast by myself are welcome, if there is no available solution.

kind regards

I asked Google again for you, but I did manage to find some hits in the end. There is already a very good question on stackoverflow concerning the exact same thing you are interested in.

There is a very nice explanation of split-and-merge provided in one of the answers, as well as simplified pseudocode.

The other answer provides a link to the implementation (although I did not check it, I'm not sure which language the implementation is in).

If you want to use split-and-merge as a reference segmentation technique to compare your idea with, other segmentation methods offered by OpenCV might also be interesting:

functions useful in image segmentation from the OpenCV manual
a chapter dealing with image segmentation in OpenCV from Gary Bradski; Adrian Kaehler: Learning OpenCV
article about Pyramid segmentation
Mean-shift segmentation: stackoverflow answer, scientific article
Watershed segmentation seems fairly well explained by Wikipedia

In the end, you can choose one of the already offered methods as the reference method after all, or, if you want to present amazing, all-encompassing comparisson, you might also want to use more than one reference method


Classical spectrum analyzer model I am trying to reproduce a numerical model of classical spectrum analyzer, described in "Fundamentals of Spectrum Analysis" by Christoph Rauscher.
Classical model consists of a few local oscillators (by number of frequency sweep points). The LO signals should be mixed with analyzed signal one-by-one, and then products should be convoluted with low pass filter impulse response. In the end, power of each filtered signal should be calculated, thus one can plot these values vs. swept frequencies.
In spectrum analysis the most popular IF band pass filter is Gaussian, so i use Gaussian low pass as equivalent.

The resulting spectrum looks bad (see figure). 3 dB bandwidth is not reproduced. Gaussian form distorted near carrier frequency.
Any suggestions?
Matlab MWE:
clear variables;
close all;

%% Signal properties
fc = 1e6;                               % Carrier frequency, Hz
fs = 10e6;                              % Sampling frequency, Hz

%% Spectrum analyzer specification
span = 100e3;                           % Span window, Hz
center = 1e6;                           % Center frequency, Hz
sweep = 1e-3;                           % Sweep time, s
rbw = 10e3;                             % Resolution filter bandwidth, Hz
points = 501;                           % Sweep points

%% Signal synthesis
samples = sweep*fs+1;

% Time scale
time = linspace(0, sweep, samples);

% The signal
signal = cos(2*pi*fc*time);

%% Gaussian filter design
order = ceil(fs/2/rbw)*10;
alpha = 10;
cutoff = rbw/2;
window = gausswin(order+1, alpha);
b = fir1(order, cutoff/(fs/2), 'low', window, 'scale');

%% Signal analysis
start = center - span/2;
stop = center + span/2;
freqs = linspace(start, stop, points);
power = 10*log10(sum(signal.^2)/samples*(samples-1-order))-3;

product = zeros(samples, points);
values = zeros(1,points);
for i = 1:points
    product(:,i) = (signal.*cos(2*pi*freqs(i)*time))';
    product(:,i) = filter(b,1,product(:,i));
    values(i) = 10*log10(sum(product((order+1):end,i).^2))-power;
end;

plot(freqs,values);
xlim([start stop]);
ylim([-120 0]);
grid on;

I believe the issue is the OP is simulating a spectrum analyzer with a "Zero-IF" architecture so is seeing the image issues as the sweep signal approach the input signal (since both are real signals). The product is down-converted to baseband with a real LO and passed through a low-pass Gaussian filter, so the resulting signal as it approach the DC center as a very low frequency real sine-wave has increasingly fewer cycles (and ultimately partial cycles) which then contributes to high variability in the power estimate. By using a complex LO instead, the envelope of the resulting signal is constant eliminating that variability (or by observing the two real components of the baseband complex signal $e^{j\omega t} = \cos(\omega t) + j\sin(\omega t)$ we see that as one partial cycle is decreasing the other is increasing to maintain the constant power at the detector).
I confirmed this was the case by changing the LO to a complex signal in the OP's code which would then match the traditional architecture of an actual Zero-IF receiver:
product(:,i) = (signal.*exp(-j*2*pi*freqs(i)*time))';

and use a complex conjugate product for the power (which also could be done using abs(x).^2)
values(i) = 10*log10(sum(product((order+1):end,i).*(conj(product((order+1):end,i)))))-power;

Resulting in:

The other option (that I would NOT recommend for simulation! - unless you need to prove to yourself that it is identical) is to move the filter from being a low-pass filter to a band-pass filter and offset the LO such that an IF frequency passes through the filter (matching traditional spectrum analyzer architectures).
For simulation purposes I would always simulate the equivalent baseband signal using complex signals rather than actual carriers which would be significantly more time consuming to process.
Further, this process can be done in one operation using an FFT which would be significantly faster and for a stationary input such as this produce an identical result. The Gaussian filter could be implemented by multiplying the time domain waveform by the Gaussian window prior to taking the FFT. For non-stationary signals further additions could be made to properly emulate the effect of the sweep time.
As far as the bandwidth, the Gaussian window used in the OP's code with
$\alpha = 10$ has an equivalent noise BW of 5.64 bins as computed from:
$$\text{ENBW} = N\frac{\sum (w[n]^2)}{(\sum w[n])^2} \tag{1} \label{1}$$
That said for a sampling rate of 10e6, and a target resolution BW of 10KHz, this would require a time capture of
$$N = 5.64 \times (10e6/10e3) = 5640 \text{ samples}$$
This would then also match the duration for the Gaussian filter for a 10KHz resolution BW in the OP's code by selecting the coefficients directly from the window function of that length, as in:
b =  gaussian(5640, 10)

I use equivalent noise bandwidth (ENBW) as resolution BW instead of the 3 dB BW, as it can be directly computed from the window as in \ref{1} and accurately indicates the same power that would be in a brick-wall filter of that bandwidth for a white noise source. The ENBW is slightly higher than the 3 dB bandwidth (on the order of 5% but depends on the window).
As recently reiterated by fred harris at the DSP Online Conference (https://www.dsponlineconference.com/), the Gaussian window is NOT the best choice. As he detailed, the reason Gaussian is often selected for spectral analysis is that it is believed to have the minimum time-bandwidth product, so minimum localization in time and frequency. However this is only true for a true-Gaussian which requires an infinite time support and only in that case when $t \rightarrow \pm\infty$ results in the minimum achievable time-bandwidth product BT of $BT = 0.5$. When we truncate the Gaussian response for the time-limited window (as required for practical implementation) the truncated Gaussian $BT>0.5$ and is no longer the minimum BT window function available for the case of finite time duration.
A window that has minimum time BW product under finite time support is the DPSS (Digital Prolate-Spheroid Sequence) window and very close to that is the Kaiser Window, both of which are available in MATLAB, Octave and Python scipy.signal. These windows also have decreasing stop-band attenuation which is also attractive for many applications. For more info on this see: https://ccrma.stanford.edu/~jos/sasp/Slepian_DPSS_Window.html , https://www.dsprelated.com/freebooks/sasp/Kaiser_DPSS_Windows_Compared.html#fig:dpsstest and https://ieeexplore.ieee.org/document/543677


What is Wavelet packet transform (WPT) time complexity I have read in many papers that time complexity of WPT is $O(N \log(N))$, but I would like to know what is the effect of Levels in the time complexity.

Must it be one of the following ?

$O(2^{L} ~N~ \log(N))$
$O(2^{L+1} ~N~ \log(N))$
$O(L ~N ~\log(N))$

$ L$ is the number of decomposition levels

Can anyone helps which is true?

I just found the answer in a paper. so thought it might be useful for someone.
If we decompose a signal with length of $N$ to $L$ level then time complexity is
$O(L\times Nlog_{2}N)$


Calculating cross-correlation using Walsh-Hadamard transform I am trying to implement MLS method of measuring impulse responses. There is an article describing the method: http://www.commsp.ee.ic.ac.uk/~mrt102/projects/mls.... As I understand, to get an impulse response of a system I need to calculate cross-calculation  of input signal - MLS and output signal of a system.
I tried to do measurements on simple lowpass RC-filter. The result I get is strange, it is a noisy signal that does not look like lowpass filter signal response. I tried to use convolution and circular correlation but it didn't help.
I suspect the problem is related with calculation methods. I calculate correlation using the convolution theorem that allows to calculate convolution using Fourier transforms of signals. But, the article describes Walsh-Hadamard transform to use in the method.
Also I checked the form of MLS, the generated signal is correct, but I am not sure should it be periodic or not.
My questions are:
Can I get wrong results because I don't use Walsh-Hadamard transform during calculations?
Could you, please, give a sample of code in Matlab or Python that calculate correlation using Walsh-Hadamard transform?
Can there be other sources of problems, for example, distortion related with periodicity of MLS-signal?
Thank you in advance.
P.S. I can show screnshots of wrong impulse responses I got if it is required.

The link provided as a reference is broken. But, regardless of what the OP's source says, the Walsh-Hadamard transform method (multiply the data vector with a $2^n\times 2^n$ Hadamard matrix $H_n$ in Sylvester form $$ H_n = \left[\begin{matrix} H_{n-1} & H_{n-1}\\ H_{n-1} & -H_{n-1}\end{matrix}\right]$$ does not give the cross-correlation of the data vector with a binary maximum-length linear feedback shift register sequence (presumably the MLS referred to by the OP) of period $2^n-1$. It is possible to permute the data vector  entries to get a close approximation of the result, indeed even the exact result, of the desired cross-correlation, but just applying the Walsh-Hadamard transform to the unpermuted data vector results in garbage, as the OP has discovered.
A description of how the MLS method works can be found in this answer of mine. The necessary permutations to make the Walsh-Hadamard transform method work are probably there somewhere on this forum.


Should we use all $N$ correspondence points to compute Homography or just $4$ points? Given correspondences of $(x,x')$, we want to find a Homography matrix $H$ that maps $x' = Hx$. If we use 4 corresponding non - collinear points and the matrix $H$ has rank $n-1$, then the homogenous equation $Ah =0$ has a unique solution up to a constant. Suppose we find more than $4$ corresponding pairs of points, we get an overdetermined system of equation and we have to use least squares minimisation to find the solution.
I am curious as to whether it is better to just use $4$ points or $>4$ points to find the Homography ?
Also, it makes sense that from a linear algebra perspective, I am unable to solve the homogenous equation exactly when we have an overdetermined system. However, from a planar perspective, given my correspondence points are good, shouldn't there be an exact solution ?

It's generally better to use all the points you have, since often, errors are independent per point.

given my correspondence points are good, shouldn't there be an exact
solution

Sure, but in realistic situations, it's not usually the case and/or you don't know that the points are good / which points are good.
A typical approach is to use all the points you have, but reject "bad" points using outlier-rejectiong  (e.g. RANSAC ) These slides have an example of such a procedure.


Getting loudness of a track with RMS I'm trying to calculate the loudness of an audio track I have stored in a buffer. The buffer contains PCM data of the signal and I want to get how 'loud' it is by using Root Mean Squared. I assume I can do this in the time-domain instead of having to switch to the frequency domain. What would be the pseudo-code for doing this?

Would I simply sample for one second (audio[0] - audio[44099], audio[44099] - audio[88199] etc..) and calculate the RMS of those values? So, for example, would I do this:

$$RMS = \sqrt{\frac{\text{audio}[0]^2 + \text{audio}[1]^2 + \text{audio}[2]^2.....\text{audio}[44099]^2}{44100}}$$

for each second?

Another thing is that the RMS value is not very well correlated with perceived loudness. You might want to consider calling it level or volume instead.
There is something called equal loudness contours which quantifies how sensitive the ear is to one particular frquency compared to another frequency, see the Wikipedia article. These curves are level dependent.
For instance, the ear is very sensitive to a 1kHz tone compared to a 100Hz tone, as shown in this image (horizontal axis is frequency in Hz):

One of the relative simple things you can do is to filter your PCM data with an inverted equal loudness curve. Or you can apply the standard A weighting, see the Wikipedia Weighting Filter article. Then you can compute the RMS value of the output of the equal loudness weighted filter.


Kalman filter restricted by areal constraints I'm developing a device for tracking pets in a room with several motion detectors. I use a Kalman filter to estimate the position, which is based on the active outputs of the motion detectors. If an motion detector is active the center of the detection region is the measurement value (position). In some cases the filter estimates the position to be outside the room.

How can I restrict the Kalman filter, should I restart and reinitialise the filter when the estimated position is not in the predefined area (of the room)?

The simplest solution in terms of computational resources as much as simplicity, I guess, is the Unscented Kalman Filter (UKF). There, the clipping might easily be added to the calculation of the sigma-points, in a way that they remain in certain boundary. I recommend that you take a look at this article, for example (working in 2016.04.22, Constrained State Estimation Using the Unscented Kalman Filter, Rambabu Kandepu, Lars Imsland and Bjarne A. Foss)


Image processing: Image recognition and classification on fish I need to identify the type of fish caught from fish images. How can I locate anchor points/landmark points to extract features from the image?

Namely, I  want to locate eye position, dorsal and pelvic fin. Need to get Fish mouth length, Dorsal and Caudal fin length.

Right now I am trying with SIFT method to get key points.

Can someone suggest me how can I get the specific key points?

Example images can be found in Fish Recognition Based on Robust Features Extraction from Size and Shape Measurements Using Neural Network, 2010, Journal of Computer Science (online version).

How can I get the specific key points?

I'm not sure but I think you are talking about Feature matching using SIFT.

Extract keypoints using SIFT from reference images of what you are trying to match.
Extract keypoints using SIFT from the test images.
Try to match the keypoints from reference images to test images using algorithms like RANSAC to get the position of the features(in your case eye,dorsal and pelvic fin).

You can find more information here.


How to calculate arbitrary phase shift in discrete signal? When I have digital audio signal, represented in samples, I can easily phase-shift the signal by 180° by simply inverting the value of each sample.

But what if I want to phase shift an arbitrary voice audio signal by say +90° or -90°?

How do I transform an arbitrary, discrete signal by phase-shifting it by an arbitrary phase angle?

The answers to this question explain the solution for the continuous-time case. In discrete time, the solution is completely analogous. The impulse response of a phase shifter with phase shift $\theta$ is given by

$$h[n]=\cos(\theta)\delta[n]+\sin(\theta)g[n]\tag{1}$$

where $g[n]$ is the impulse response of an (ideal) discrete-time Hilbert transformer. Eq. $(1)$ shows that the solution is just a linear combination of the identity system and a Hilbert transformer (which performs a $90$ degree phase shift).


How to apply a 2D -2D homography matrix to an image? I have a mini project to do about homography, but I am really a beginner in Matlab. I have calculated the homography matrix easily by taking a set pixels from both images but the problem is how to apply this matrix to all pixels of the original image (I think it's difficult to treat each pixel individually).

Here is the original image :

I want to apply the homography to this image (I have already extracted the homography matrix) to get a frontal view of the cover.

So with my weak knowledge of Matlab, I applied the following operations:

I get 4 points from the corners of the original image and the corresponding points in the output image
I have estimated the homography :

 0.9638   -0.0960   52.5754
 0.2449    1.3808  -17.0081
-0.0001    0.0013    1.0000

then I have applied it on the image using this poor Matlab code:

function [Ired] = affect(I,P)
I1=im2double(I);
for i=1:308
i
for j=1:288
v1 = [j;i;1];
v2 = P*v1;
v3 =v2/v2(3,1);
if (v3(1,1) <0 || v3(2,1)<0 || v3(3,1)<0)
   continue;
else
Ired(round(1+v3(2,1)),round(1+v3(1,1)))=I1(i,j);
end
end

end

and the result is the following:

but the problem with this code is it takes a lot of time (sometimes 10 minutes) to perform the operation on the image, especially if the size of the image is high. Is there a way to implement this code without loops to make it faster?
I get these black dots all over the image, I can do a dilation but the image will be distorted, do you have an idea how to eliminate this dots without distorting the image?

I cannot use predefined Matlab functions in this project.

What do you mean by "apply this matrix to all pixels of the image without treating them one by one"?

The image is made of pixels, each one with its x-y location. The homography maps a coordinate on a plane into another one so you actually HAVE to apply it to every single pixel

(x' ; y' ; 1) = H (x ; y ; 1)

where H is your 3x3 homography matrix. Doing this you'll get the new image (plus some black areas), simply iterate the equation for every pixel in the image.

(sorry for the bad formatting of the formula, wrote it down fast)


Calculate the magnitude and phase of a signal at a particular frequency in python I have a signal for which I need to calculate the magnitude and phase at 200 Hz frequency only. I would like to use Fourier transform for it. I am very new to signal processing. And this is my first time using a Fourier transform.
I found that I can use the scipy.fftpack.fft to calculate the FFT of the signal. Then use numpy.mag and numpyh.phase to calculate the magnitude and phases of the entire signal. But I would like to get the magnitude and phase value of the signal corresponding to 200 Hz frequency only. How can I do this using Python?
So far I have done.
from scipy.fftpack import fft
import numpy as np

fft_data = fft(signal)
magnitude = np.mag(fft_data)
phase = np.phase(fft_data)

You can find the index of the desired (or the closest one) frequency in the array of resulting frequency bins using np.fft.fftfreq function, then use np.abs and np.angle functions to get the magnitude and phase.
Here is an example using fft.fft function from numpy library for a synthetic signal.
import numpy as np
import matplotlib.pyplot as plt

# Number of sample points
N = 1000

# Sample spacing
T = 1.0 / 800.0     # f = 800 Hz

# Create a signal
x = np.linspace(0.0, N*T, N)
t0 = np.pi/6   # non-zero phase of the second sine
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(200.0 * 2.0*np.pi*x + t0)
yf = np.fft.fft(y) # to normalize use norm='ortho' as an additional argument

# Where is a 200 Hz frequency in the results?
freq = np.fft.fftfreq(x.size, d=T)
index, = np.where(np.isclose(freq, 200, atol=1/(T*N)))

# Get magnitude and phase
magnitude = np.abs(yf[index[0]])
phase = np.angle(yf[index[0]])
print("Magnitude:", magnitude, ", phase:", phase)

# Plot a spectrum
plt.plot(freq[0:N//2], 2/N*np.abs(yf[0:N//2]), label='amplitude spectrum')   # in a conventional form
plt.plot(freq[0:N//2], np.angle(yf[0:N//2]), label='phase spectrum')
plt.legend()
plt.grid()
plt.show()

And here is a useful manual with detailed explanations: reference.


How much deviation is acceptable in FFT calculations (using computer algorithms)? It seems apparent that different FFT algorithms (e.g. Numpy.fft vs FFTS) produce slightly different results. I'm seeing deviations around 0.005 to 0.01 between programs.

The reason for these, I guess, lies in computer arithmetic, but are there simple ways for finding out where the differences are caused? Or to tell, which algorithm is more precise?

I wouldn't like to give a trivial answer as "it depends on the application" but indeed it does depend on the application!

Considering the fact that FFT computation will not be the only source of error in many applications it is reasonable to allow some acceptable deviations but generally being bounded by such as the required precision of the overall system.

On the other hand very strict precision bounds may be defined on such as IEEE floating point hardware (or soft simulation) units so that a prescribed set of computations will not generate an error larger than very small epsilon values.

I don't know if such a set of normalized limits exist specificaly for FFT or like algorithms. If you find out that there are mismatches from ideal DFT values and the corresponding fast algorithms you must make your own decision as to whether you should allow this error in favor of some other gains.


How can we measure the quality of the FIR filter? The process of designing filters includes too many parameters like bandwidth, transition, degree of attenuation in side-bands, phase delay, causality, the magnitude of ripples, and so on. First of all, we want to approach the ideal filter structure so we can design many different filters to obtain the desired ideal filter.
How can we compare all of these filters in terms of quality? Is there any way to measure the performance of filters? I am looking for a mathematical solution to this problem.

The answer can be given without being at all specific to filters. In fact, one has to answer that in terms of engineering in general:
Whenever you have alternative solutions, the only measure of quality that really matters, is fulfillment of the requirements.
So, you'll have to know what you'll use that filter for, and evaluate how well the alternatives perform with respect to that specific application. Nothing else helps you. There's arbitrarily many properties that a filter can have, and for most of them, it's not clear whether a smaller or bigger value is better, and they all are intrinsically linked to other parameters. There's no single "best" solution, all depends on the application.
So, define your application. Define how that sees performance. For example, if your FIR is used in a wireless communications system, then the FIR that makes the system have the least bit error rate might be the "highest quality" one for this one system under this specific metric that you care about. For a different system, a different FIR will be the best one.


Impulse Response with non-delayed feedback I've been trying to solve a Discrete LTI System problem where the overall impulse response of the following System, whose input and output sequences are, respectively, $x(n)$ and $y(n)$, has to be found:

This is what I've tried:

Calling the output sequence that comes out of the adder $e(n)$, I know that $y(n) = e(n) \star h(n)$. I also know that $e(n) = x(n) - y(n)$, which leads to $y(n) = x(n) \star h(n) - y(n) \star h(n)$.
Thats the point where I got stuck. I don't know how to proceed in order to isolate $y(n)$ in the last expression. I would like to isolate it in order to make $x(n) = \delta(n)$ and evaluate the overall Impulse Response.

Also, I am not supposed to use any transform techniques as the book I am following has not yet arrived at those subjects. If possible, don't show the final calculations, only what I should do in order to separate $y(n)$. I would like to figure the rest out by myself.

Thanks in advance.

Take the z-transform of your equation to get $y(z) = x(z) h(z) - y(z) h(z)$ which you can simplify to $$ y(z) \left( 1 + h(z) \right) = x(z) h(z) $$ and therefore get the transfer function
$$ H(z) := \frac{y(z)}{x(z)} = \frac{h(z)}{1+h(z)} $$

Now plug in the z-transform $h(z)$ of your response $h(n)$ and perform an inverse z-transform to get the total impulse response $H(n)$.

An approach that is more basic and does not require z-transforms works for your particular example by solving the recursion equation explicitly. The inner system with the response $h(n) = u(n) \left( \frac{1}{2}\right)^n$ can be written with the recursive difference equation $y'(n) = x'(n) + \frac{1}{2} y'(n-1)$ where $x'$ and $y'$ are the output and input of the subsystem. Then we can write the full system in terms of the following equations:

$$ y(n) = y'(n)$$
$$ x'(n) = x(n) - y(n)$$
$$ y'(n) = x'(n) + \frac{1}{2} y'(n-1) $$

Eliminating $x'$ and $y'$ results in

$$ y(n) = x(n) - y(n) + \frac{1}{2} y(n-1) $$

which simplifies to

$$y(n)=\frac{1}{2}x(n)+\frac{1}{4}y(n-1)$$

This system has the impulse response $$H(n) = \frac{1}{2} u(n) \left(\frac{1}{4}\right)^n$$


Recovering a signal that was upsampled by non-integer factor I have a data sequence $G[n]$ that is upsampled at the transmitter of a communication system by a non-integer factor ($1.5$ to be precise). At the receiver, I am supposed to downsample the data and recover the sequence $G[n]$.

In MATLAB, I am implementing both the transmit and receive signals. When I upsample the signal $G[n]$ by integer factors at the transmitter and then downsample the received signal by the same integer factors, I am able to successfully recover the exact $G[n]$, and continue further processing.

However, when I upsample $G[n]$ by a non-integer factor $1.5$, and then downsample the received data by the same factor $1.5$, I do not get the exact sequence $G[n]$ that was transmitted. In MATLAB, for the transmitted sequence, I simply used commands upsample with factor $3$ and downsample with factor $2$. At the receiver I upsampled by factor $2$ and then downsampled by factor $3$.

Is there a way to recover the exact sequence $G[n]$ at the receiver even if it was upsampled by a non-integer factor at the transmitter?

~ryan

The Sampling operation (both upsample and downsample) depends on two very critical conditions:

1- The existance and applicability of ideal frequency selective filters

2- The operated signal being strictly bandlimited

For most practical systems these two critical conditions are only approximately met. Hence the computational results by those practical systems are only approximate and not exact.

Note that in case of simple but crude interpolation (upsampling) by zero-order hold, which requires no filter at all as being just a sample replication, you can exactly recover the original samples due to the simplicity of the operation.

Unless some advanced algorithms are utilized, therefore, some deviation (but controllable) is expected at the output of any practical sample rate conversion systems.