texts:
I don't understand the meaning of scaling function when I used CWT/DWT.
For me, I understand that the wavelet function (mother wavelet) is the function that we used to decompose our signal in subsignal corresponding to the frequency of the mother wavelet (for example db8). (Maybe DWT is different from wavelet decomposition ?)
For example, if I decompose my signal into 5 levels (with wavedec in MATLAB), I have :
D1, D2, D3, D4, D5 and A5, with D the detail's signal and A the approximation's signal.
But here I never used the scaling function to make the decomposition of my signal.
Can you help me on how to understand the scaling function ?


texts:
Assuming normalized biquad coefficients:

$$
H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}} {1 + a_1z^{-1} + a_2z^{-2}}
$$

And taking the definition of flipping the polarity as multiplying every tap in a signal by -1 i.e.

$$[0.2, 0.1, -0.1, -0.2] \to [-0.2, -0.1, 0.1, 0.2]$$

For a zero delay biquad is it sufficient to check that $b_0 >= 0$ to assume the biquad is not flipping the polarity accidentally (although it will be modifying the waveform) and if it were, a correction $H'(z)$ would be:

$$
H'(z) = \frac{-b_0 - b_1z^{-1} - b_2z^{-2}} {1 + a_1z^{-1} + a_2z^{-2}}
$$

This relates to audio.


texts: Output gain when upsampling and downsampling

texts:
Do the processes of upsampling and downsampling affect the magnitude of the transform of a signal? And if not, why am I seeing everywhere that a filter with gain different from 1 is applied after up/downsampling a signal?

EDIT: For example in Discrete-Time Signal Processing by Oppenheim-Schafer (Section 4.6, 2nd edition):

Also, in the summaries of the subject that my teacher did himself, up- and downsampling are always followed by a lowpassfilter with gain different from 1.


texts: How to decide the number of cepstral coefficients?

texts: For audio classification we are using cepstral coefficients, but how should I decide the number of coefficients?


texts: Signal reconstructed using IFFT: amplitude scale larger than input

texts:
I am trying to test reconstruction of signal from IFFT.
Below is the python program. I break the signal into overlapping windows and do FFT after multiplication with window function. Then I apply IFFT to the FFT output and then overlap-add to reconstruct the original signal. My problem is the reconstructed signal's amplitude is almost 10 times the input. Why is this so ? How can I scale it ?

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

winsize=512
shift=16

aircutrawdata=pd.read_csv('input.csv',sep=',',skiprows=8000,skipfooter=3500).values
#%%
vibaxisdata=aircutrawdata[:,8]
data_len=len(vibaxisdata)
window_func=np.hanning(winsize)
fftdata=np.zeros((0,int(winsize/2)))
rec_frame=np.zeros(len(aircutrawdata))
for frame in range(0, data_len, shift):
# =============================================================================
#     if frame>0:
#          break
# =============================================================================
    kiri=vibaxisdata[frame:frame+winsize]
    if len(kiri) != winsize:
        break
    windata = window_func * kiri
    fftframe=np.fft.fft(windata,n=winsize)
    ifft_frame=np.real(np.fft.ifft(fftframe,n=winsize))
    rec_frame[frame:frame+winsize]+=ifft_frame

fig,(ax1,ax2)=plt.subplots(nrows=2,ncols=1)
ax1.plot(vibaxisdata,color='b')
ax2.plot(rec_frame,color='r')
plt.show()

The wave on the top is the original input and the wave at the bottom is the reconstructed signal. You can see that the signal is rather faithfully reconstructed except the amplitude. It is almost ten times the original. Is that how it is supposed to be ?


texts: What does the uncertainty principle say about recursive filters?

texts:
The uncertainty principle is usually stated as a relationship between a continuous signal and that signal's Fourier transform, and says that
$$
\int_{-\infty}^{\infty} \! x^2 f(x) \ \mathrm{dx} \int_{-\infty}^{\infty} \! \xi^2 \hat{f}(\xi) \ \mathrm{d\xi} \ge \frac{1}{16\pi^2}.
$$

Apparently it is somewhat difficult to define an equivalent relationship for discrete time signals, but, for "reasonable" bandlimited signals with "reasonable" transforms, it is basically correct to say that the above holds true of the continuous signal reconstructed from the samples.  For a good overview I found the following non-paywalled paper useful:
Venkatesh, Kumar Raja, Vidyasagar: On the uncertainty inequality as applied to discrete signals, Int'l J. Math. and Math. Sci., 2006:48185, 2006. doi:10.1155/IJMMS/2006/48185.

So I also understand what this says about FIR filters: I can treat my kernel as a signal, apply the uncertainty principle and see that if I want a very narrowband filter, I'm going to have to have a lot of taps.

My question is what does the uncertainty principle say about an IIR filter?  If I have a very narrow band IIR filter it would seem that the uncertainty principle would say that I would have to trade something off.  What?  I can't figure out what the appropriate signal is that I should take the variance of.  Is it continuous reconstruction of the impulse response of my filter?  (I.e., so I should expect a very narrow-band IIR filter to have an extremely long startup transient?)


texts: Determining the mean and standard deviation in real time

texts:
What would be the ideal way to find the mean and standard deviation of a signal for a real time application.  I'd like to be able to trigger a controller when a signal was more than 3 standard deviation off of the mean for a certain amount of time.

I'm assuming a dedicated DSP would do this pretty readily, but is there any "shortcut" that may not require something so complicated?


texts: Correct the phase of several damped sines

texts:
I have a signal $s(t)$ which consists of 3 damped sines it can be written as:
$$
s(t) = A_1\sin(\omega_1t + \phi_1)\exp(-t/t_1) + A_2\sin(\omega_2 t + \phi_2)\exp(-t/t_2) + \\A_3\sin(\omega_3 t + \phi_3)\exp(-t/t_3) + n(t),
$$
where $n(t)$ is noise.

How to make the $\phi_1 = \phi_2 = \phi_3 = 0$; Any filter can fulfill the task?
It is also desirable to extract $A_1$, $A_2$ and $A_3$. But some times $A$s can be negative.


texts:
I have seen several methods to calculate Autocorrelations using FFTs, and am confused about why they differ.

Zero-Pad it to double its original length.Take the FFT. Then replace all the coefficients with their square magnitude, except for the first, which you set to 0. Now take the IFFT. Divide every element by the first one. This gives you the autocorrelation. https://stackoverflow.com/questions/4583950/cepstral-analysis-for-pitch-detection

create a centered version of x by setting x_cent = x / mean(x);
pad x_cent at the end with entries of value 0 to get a new vector of length L = 2^ceil(log2(N));
run x_pad through a forward discrete fast fourier transform to get an L-vector z of complex values;
replace the entries in z with their norms (the norm of a complex number is the real number resulting of summing the squared real component and squared imaginary component).
run z through the inverse discrete FFT to produce an L-vector acov of (unadjusted) autocovariances;
trim acov to size N;
create a L-vector named mask consisting of N entries with value 1 followed by L-N entries with value 0;
compute the forward FFT of mask and put the result in the L-vector adj
to get adjusted autocovariance estimates, divide each entry acov[n] by norm(adj[n]), where norm is the complex norm defined above; and
to get autocorrelations, set acorr[n] = acov[n] / acov[0] (acov[0], the autocovariance at lag 0, is just the variance). https://lingpipe-blog.com/2012/06/08/autocorrelation-fft-kiss-eigen/

And of course there is simply the first method but without setting the first element in the frequency domain to zero (i.e. what you'd expect from the Wiener–Khinchin theorem) https://stackoverflow.com/questions/3949324/calculate-autocorrelation-using-fft-in-matlab.
Why do the first two methods differ from the naive implementation of Wiener–Khinchin theorem. Are there any references that may help me figure out which method is best (and maybe on fast convolutions in general)?


texts: FFT time domain average vs frequency bin average

texts:
I have multiple trials of physiological data. I am doing a frequency based analysis to analyze power (amplitude) in certain frequencies of interest. Is averaging multiple trials of equal length and then taking a single FFT of the averaged signal vs computing FFT for each trial and then averaging the frequency bins the same? In practice im finding this not to be the case.

Specifically, the signal naturally has a strong 1/f component and this gets emphasized if I compute the FFT of each individual trial and then average the amplitudes (real part) of each frequency bin. Are the two equivalent? is there a right way to do things? or under what principled conditions should the choice between time domain averaging vs frequency bin averaging be made?


texts: Introduction to Clustering and Classification of Linear Sequence

texts:
I'm a programmer/engineer that is getting into developing prototype wearable technology. I would like to be in a position to interpret raw accelerometer data and categorised the data into activities like running/walking/rest, etc.

The building and coding of the hardware and software is relatively straight forward, but I do not know the terminology and mathematical processes required in order to take the raw data and turn it into a training set that can be used for a real time system.

My question is this: What is a good DSP starting point to take a time based linear series and extract features so that they may be identified at a later stage?

At present, I have an tri-axis accelerometer device which samples at 50Hz with a range of -/+4g. I would like the capability of being able to readily identify activities based upon a training set, or machine learning process, if that's even possible?

Apologies that this isn't necessarily a specific DSP question, thanks for you assistance.


texts: 'imread' function in MATLAB gives image data.what this image data is?is that intensity of pixels?and how MATLAB computes this data?

texts: i checked imread function code but couldnt understand how it is calculating the data from given image.Please explain.


texts: All-pass phase compensation only in a specific frequency range

texts:
I am trying to design an all-pass filter in order to compensate the phase distortion at the output of an equalizer. To do so, I take the phase vector of the output, compute the real and imaginary part using  $e^{j\phi}$, then I do an IFFT and save the time-reversed output. The saved impulse response of the all-pass filter looks like this:

Problems start when I try to restrict the working range of the all-pass equalizer: for example, let's say that the phase compensation needs to work only in 400-16k Hz range, leaving the remaining frequencies uncompensated. To do so, I'm essentially zeroing the phase outside the working range before computing the real and imaginary parts. Probably this is not the correct way, since the resulting impulse response has several ringings at the start and the end of the plot, as you can see from the following picture.

Could you please help me to understand which is the problem here?
EDIT in response to @Dan Boschen's answer
This is the solution I implemented in Octave following your suggestions:
h = read_dat("heqm.dat"); % the IR to equalize
NFFT = length(h);

H = fft(h,NFFT);
ph = unwrap(angle(H));

x = 1:NFFT;
p = polyfit(x,ph,1);
f = polyval(p,x);
f = f-f(1);
figure; plot(ph); hold on; plot(f);

The corresponding plot is:

EDIT 2
The real(ifft(exp(1j*ph))) is:

Next, the definition of the mask:
f1 = 100;
f2 = 300;
half_mask = zeros(1,NFFT/2+1);
half_mask(f1:f2) = 1;
mask = [half_mask, fliplr(half_mask(2:end-1))];
idx = find(mask);

res = f;
res(idx) = ph(idx);
figure; plot(res);

which corresponds to:

Finally, the IR computation:
imp = real(ifft(exp(1j*res)));
figure; plot(imp)

which unfortunately has aliasing:

Did I do something wrong?


texts: Question regarding DTFT of a complex signal

texts:
I have been doing DTFT practice problems for my DSP course, and I encountered this problem in the textbook that completely stumped me.

The question asks to find the DTFT of the shown signal and to plot its magnitude and phase. I am very comfortable with doing this with real-valued signals, but have no idea where to start when it comes to complex signals. Especially how the plotting would be done. I would appreciate any advice on how to tackle the problem.


texts: Differentiating between zero and not sending for OOK

texts: How would a trigger works for on-off keying (OOK)? How would it differentiate between not receiving anything and zero?


texts: Zero padding effect on a FFT of gaussian noise

texts:
I have a gaussian noise $\nu(t)$ with variance $\sigma^2$. After a FFT I get $X(\omega)$. If now I do the IFFT on the $X(\omega)$ can I say that the result is still a gaussian noise of variance $\sigma^2$? What is the effect of zero padding on FFT? How the zero padding affects the statistics of the noise after FFT and IFFT?
Thanks.


texts: Calculate ifft using only REAL forward FFT

texts: There are already questions on this topic here and here. I am trying to calculate the ifft using a DSP that only has a library function for a forward REAL (assumes conjugate symmetry) fft. For a given a signal of length N, the output of the FFT process is length N interleaved real/complex values, where the real part of the nyquist bin is put into the imaginary part of the DC component (always zero for a real signal - so this keeps a constant memory footprint). I'm having trouble fathoming how to calculate the inverse using such a function.


texts: Bandpass Filter window length: filter "in phase" with sinusoidal input

texts:
Say you have a bandpass filter and you place it on 40 discretely sampled points from a sinusoid with a period of 40.

Why is the bandpass filter (green curve on the image) only "in phase" with the sinusoid if you use a window length (impulse response length) for the filter of 1/4th the period of the sinusoid (here 40/4=10)? So the green filter has the same turning and bending points only if you use the 10 past observations, multiply with the impulse response (of length 10) and add that up..

Am I correct that this is related to how the length of the impulse response affects the phase delay?


texts: Optimal pulse shape for minimal interference in adjacent frequencies

texts:
I need to find an optimal pulse shape, that convoluted with matched filter on adjacent frequency will give minimal convolution/correlation (Minimal interference) for a given pulse length.

For example Pulse Length of 10 msec, 200 Hz separation (10,000 Hz and 10,200 Hz)

Image below demonstrating the problem.  Currently I have interference with adjacent signal of approx 0.05.

What would be the Optimal pulse shape for given pulse length?


texts: Compute the time derivative of a noisy digital signal?

texts:
The issue is that my signal is very noisy. I need extract its time derivative as accurate as possible.  P.S. I do not have any prior knowledge on the signal (black box).
On forums some suggested Savitzky-Golay filter.
Any idea please? If so, is there any c++ library for the purpose?
In fact, for my application I need to compute optical-flow like information for control purpose. I compute an estimate using image information. Then I need to compute the time derivative of this estimate.
4th order Savitzky Golay filter introduces delay, yet I need the output in real-time (real time control). For info:

The signal is regularly sampled;
The noise is not defined but bounded;
The output needs to be real-time: delay-minimal;
The signal is black box: I only get a measure each iteration.


texts: State of the art algorithms for estimating noise variance

texts:
Suppose I have a signal $X[n]= f[n] + W[n]$, where $W$ is uncorrelated Gaussian white noise satisfying $E\{W[n]W[k]\} = \sigma^2 \delta[n-k]$.

What are the state-of-the-art algorithms for determining the noise variance $\sigma^2$ from $X$? (i.e., without oracle information?)


texts: Second Order Transfer Function : Imaginary component in $R$?

texts:
I have been looking at the following transfer function:

$$ H(z) = \tfrac12 -  \tfrac12 z^{-2} $$

Given the usual method for finding $\theta$ and $R$ in the complex plane, I calculate that $\theta = \frac{\pi}{2}$ and $R = j\sqrt{\frac12}$

This would initially make me think that the conjugate zero pair here is 90° from the real axis, however this imaginary $R$ component has thrown me completely - the equivalent frequency response also seems odd to me as it shows zero at DC and Nyquist with a gentle curve upwards in between.

What is going on here? I think that $R$ in the complex plane can't be imaginary, but don't understand how to arrive at the correct angle and magnitude for the conjugate poles.


texts: New to DSP and having trouble filtering data

texts:
I am trying to make program in matlab that counts the number of times the accelerometer was lifted up and down (a rep counter for exercise). Basically my phone is getting accelerometer data on all three axes and then i import this data into matlab. I have tried to filter the data using a simple filter a = [1 -0.9] b=0.1 filter(b,a,data) but that doesn't really get rid of all the high frequency noise.

Another problem i am running into is that when i apply certain filters to the data the data gets shifted so i can no longer integrate it to get velocity (my original plan was to simple filter out high frequencies and then integrate the data to get velocity and then count the number of times it crosses zero and divide by two (which would correspond to the number of times it was at a peak and trough - divided by two for number of cycles).

But the data is so noisy that when i apply a filter not only does it get shifted away from the x axis making it impossible to integrate it, it also still has noise which would give me incorrect zero-crossings.

Is there a good beginner introduction to signal processing and accelerometer data processing? What techniques should i be aware of?

If you want to try to give the data a go here is some sample data, the rep count should be 5, The data was all taken at a 50Hz sample rate. I have to put data in the comments because SE doesnt allow posting multiple links a new poster


texts: proving that a log gabor filter has 0 DC offset

texts: I have read somewhere online that the log gabor filter has an advantage over the gabor filter, in the sense that it has 0 DC component. How do you prove this property mathematically? Thanks in advance.


texts: In What Way Is the Difference of Gaussian (DoG) More Tunable than the Laplacian of Gaussian (LoG)?

texts:
I was reading through this post "What is the difference between Difference of Gaussian, Laplace of Gaussian, and Mexican Hat wavelet?" and noted in Chris Luengo's answer that he mentioned that the Difference-of-Gaussians (DoG) had an additional property of being tunable in specific way different to that of the Laplacian of Gaussian (LoG) - I've been a lurker for quite some time, so I don't have the privileges to comment, unfortunately.

My intuition is that the DoG has two sigma parameters (standard deviations) that can be varied independetly, wheres the LoG only has one (the initial Gaussian), and that this gives the DoG some bandpass flexibility that the LoG does not have. Have anyone looked into this further? Is there any research dealing directly with this issue? What I want to know are some more specific details in which way the DoG is tunable compared to the tuning properties of the LoG.

Thank you for your attention.


texts: What is the suitable mathematical form of the multipath channel impulse response?

texts:
I am trying to understand the nature of a wireless multipath channel impulse response but I am a bit confused
while interpreting the time domain multipath impulse response.
There are different mathematical models of the wireless multipath channels in some of the papers.

As far as I know the mathematical form of the time domain multipath impulse response is modelled in the following equation as:

\begin{align}
h(\tau, t) &= \sum_i a_i(t)\delta(t - \tau_i (t))\tag{Eq. 1.1}\\
H(f; t) &= \sum_i a_i(t)\displaystyle e^{-j2\pi f\tau_i (t)}\tag{Eq. 1.2}
\end{align}

Here is the link for this channel model.

So in this model; $N$ is the number of paths, $\tau$ is the delay corresponding to each path,and there is no doppler introduced in the channel. Only the path amplitudes and the path delays are changing for each symbol transmission. Amplitudes and delays are assumed to be constant during one symbol transmission.

If we want to get the frequency represantation of this channel then it takes the form in frequency domain in ($\rm Eq.\ 1.2$). Then we obtain the corresponding frequency response of this channel by giving values within the frequency range which we are dealing with.

But in the equation ($\rm Eq.\ 2.1$) below the same multipath channel is modelled in the time domain as:

\begin{align}
h(t)&=\sum_{k=1}^K a_k\cdot e^{j\phi_k}\cdot e^{-j2\pi f_c\tau_k}\delta(t - \tau_k)\tag{Eq. 2.1}\\
H(f)&=\sum_{k=1}^K a_k\cdot e^{j\phi_k}\cdot e^{-j2\pi f_c\tau_k}\cdot e^{-j2\pi f\tau_k}\tag{Eq. 2.2}
\end{align}

Here is the link for this model.

This makes me confused because this second model has the phase shift $2\pi f_c \tau$ and I don't understand why this term has came into place here.
Because this term, $2\pi f_c\tau$, appears in the form of the frequency response of the channel shown in the first equation,($\rm Eq.\ 1.2$).
But here it is a parameter of the time domain represantation.

Additionally the second model has the phase shift $\phi_k$ which I think this is a phase shift caused by some of the scatterers in the channel environment.

Could you please help me answering the following questions:

What is the difference between the time domain impulse responses of the channel models shown in ($Eq.\ 1.1$) and ($Eq.\ 2.1$)?
In order to get the time domain represantation of the channel paths (as $a+jb$), do I first get the frequency response of the channel within a specific frequency range which I want to transmit my data over those frequencies?

For example: If I want to transmit an OFDM data over 64 subcarrier frequencies, $900\ \rm{MHz}, 901\ \rm{MHz},\ldots, 963\ \rm{MHz}$) through my multipath channel, with the command filter in MATLAB, do I first need to find the channel frequency response for all 64 frequencies then its time domain values by performing IFFT, then applying the data to the filter (channel impulse response)?
timeDomainRxsignal=filter(channelImpulseResponse,1,txData)

I know there is a command in MATLAB for Rayleigh or Rician channel models but here I would like to make this channel manually for my practice.


texts: Find the valley in the given in a sequence of numbers

texts:
I am trying to find the valley in a sequence of numbers of around 2000 length.

Always prominently there will be two peaks. But the value and location might change. I need to find the valley, defined as the lowest point between the two peak values.

I tried using the following code:

% sums = list of positive integers generated from an analysis
inverted_sums = double(max(sums) - sums); % inverting the sequnce
[~,locs] = findpeaks(inverted_sums) % finding the peaks

But the above code returns many valleys with even a small difference, including the correct one. How do I solve it?


texts: Could 2 synchronised RX chains be used to receive a signal too wide for just one of them?

texts:
Say the goal is to demodulate a "wide" signal, e.g. 300 MHz wide. An SDR is available with a maximum RX bandwidth of 150 MHz (ADC sampling freq. limitation), but it comes with 2 RX channels sharing the same clock source.

The proposed approach is connecting the same antenna to both RX chains. Ideally, they would take the same clock source, but one ADC would trigger on a rising edge and the other on the falling edge of the clock, thus effectively doubling the ADC sampling rate and allowing downstream demodulation logic to correctly demodulate a 300 MHz signal.

If this is practically doable, what are the constraints (w.r.t. SDR analogue RF frontend capabilities, compensation for loss of signal strength due to splitting it into 2 RX inputs etc.) and if not, why not?


texts: What are the differences between notch filter & Goertzel?

texts:
I would like to integrate in my application a tone detector but I'm confused between 2 kinds of filters: Goertzel or Notch?
I saw that Goertzel is a recommended for tone detection. But - what if I implement a Notch filter using Matlab with order of 20 (for example)??


texts:
So, I have this problem where I want to apply a kernel to an image and count the number of matches that happened.

So for example, if I have the kernel:
$$\begin{bmatrix}
1 & 2 & 1\\
1 & 2 & 1\\
1 & 2 & 1\\
\end{bmatrix}$$
And the image:
$$\begin{bmatrix}
1 & 1 & 2 & 1 & 2\\
1 & 1 & 2 & 1 & 2\\
1 & 1 & 2 & 1 & 2\\
\end{bmatrix}$$
The result should be (ignoring the borders):
$$\begin{bmatrix}
3 & 9 & 0\\
\end{bmatrix}$$
Explanation: when the kernel is in the first position (2,2), the 3 left ones match to the kernel, in the middle position (3,2) all 9 elements match and in the right (4,2) there is no match.

I ended up doing two convolutions (or correlations), one to count the numner of ones and one to count the number of twos.

For example:
First convolution (to count where value is 1):
$$\begin{bmatrix}
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
\end{bmatrix}$$
$$\begin{bmatrix}
1 & 1 & 0 & 1 & 0\\
1 & 1 & 0 & 1 & 0\\
1 & 1 & 0 & 1 & 0\\
\end{bmatrix}$$
Result:
$$\begin{bmatrix}
3 & 6 & 0
\end{bmatrix}$$
Second convolution (to count where value is 2):
$$\begin{bmatrix}
0 & 1 & 0\\
0 & 1 & 0\\
0 & 1 & 0\\
\end{bmatrix}$$
$$\begin{bmatrix}
0 & 0 & 1 & 0 & 1\\
0 & 0 & 1 & 0 & 1\\
0 & 0 & 1 & 0 & 1\\
\end{bmatrix}$$
Result:
$$\begin{bmatrix}
0 & 3 & 0
\end{bmatrix}$$

At the end I add up the results to get the expected result. The problem is that if my image has 256 values, I have to do that 256 times and it becomes 256x slower than a convolution. Still, this ends up being faster than doing for loops and checking all kernel positions. This is partly due to the fourier transform that makes convolution faster.

My question is if it is possible to apply this operation faster? Maybe by just modifying the convolution a little bit?

EDIT:
My kernel has maximum size 64x64 and my image has size around 300x300.
If I use the standard kernel convolution (Direct Convolution), then it will try to match the kernel for every position of the image. To match the kernel it takes 64x64 comparisons, and it will do that for every position in the image (300x300), so this will take 300x300x64x64 comparisons (total 368.640.000 operations).

If I compute first the fourier transform of the images (FFT Convolution), then I can do the multiplication directly and it will cost only (300x300 operations). There is a cost to compute the fourier of the images in the first place, so at the end it will take around 300x300xlog(300)xlog(300) operations (this takes a total of around 500.000 operations).

OpenCV has a convolution method (filter2D) and a template matching method (matchTemplate) that are very fast, and I believe they use fourier plus some other optimizations.

The function matchTemplate in OpenCv has different methods to compute the match https://www.docs.opencv.org/2.4/doc/tutorials/imgproc/histograms/template_matching/template_matching.html, but none of them work for what I need.

What I need is a function to count the matches, and I want to know if it is possible to that similar to other template matching methods in opencv (that I assume use the FFT Convolution).

An approximate result should be fine, so if the output for the first example is, lets say, [3.28, 9.45, 0.32], it shouldn't be a problem.


texts: How to understand the relationship between the filter time constant and the half-power cutoff frequency

texts:
I learned that the time constant can be computed as $\frac{1}{2 \pi f_0}$, where $f_0$ is the half-power cutoff frequency of a high-pass filter.

However, I was wondering how the time constant and the half-power cutoff frequency can be related. For example, if the half-power cutoff frequency is 2 Hz, why must the time constant be $\frac{1}{2 \pi (2 \text{ Hz})}$ seconds? Why can't it be some arbitrary number, as long as at that time the output reaches $\frac{1}{e}$ of the starting value?

My knowledge of signal processing is limited to a very basic understanding of Fourier analysis and filtering, so a bit more illustration would be appreciated.


texts: Solution of Differential Equation with Complex Conjugated Variables

texts:
Consider the following first-order differential equation:
$\frac{dW(t)}{dt}+AW(t)+BW^*(t) + C=0$
where $W(t)$ is an $N \times 1$ column vector, $W^*(t)$ is its complex conjugated version, $A$ and $B$ are constant $N \times N$ matrices, and $C$ is a constant $N \times 1$ column vector.
Can anyone find the solution of $W(t)$? Thanks!


texts: Filtering out a specific Gaussian from a signal using Fourier concepts

texts:
I have been wondering about a problem and need some help.
Let's say I have a 1-d time signal (discrete), I also have the knowledge that the signal consists of repeated Gaussians with known fixed standard deviation and for starters let's say of known a fixed amplitude. So the signal consists of series of similar Gaussians and then obviously there is the random noise.
Is there any way to use that knowledge in frequency filtering to get rid of the noise?


texts:
I would like to create a set of audio samples for playback as a midi-controlled instrument. Thus: a sample for every key of the piano range, with multiple volume timbres for each note. I am able to create recordings at various specifications. My problem is finding a specification as to how a sample set is packaged so that it can be used by a standard DAW. Searching so far has not turned up any leads.

If there is a appropriate in another forum to try, I'd be happy to learn of it.


texts: Magnitude and phase of frequency response from Z-transform

texts:
I have been given a signal which I have completed the Z-transform for and have got:
$H(z)=1-z^{-1}$. I need to find the magnitude and phase of this transform but am unsure on how to actually do this. All I have got so far is $H(e^{j\omega}) = 1 - e^{-j\omega}$.


texts: The signal $\frac{e^{j10\pi n}}{10}$ is periodic/aperiodic?

texts:
How to verify the function
$$
\frac{e^{j10\pi n}}{10}
$$
is a period/aperiodic function ?

My attempt:

$\omega_o=10\pi\implies f_o=\omega_o/2\pi=5$, which is a rational number. So the given function is periodic.
The fundamental period, $N=2\pi\frac{m}{\omega_o}=2\pi\frac{5}{10\pi}=1$

Is it the proper way to approach the problem?


texts:
I'm working on an alarm recognition program in MATLAB in which multiple reference audio samples of different alarms would be recorded. Afterwards, the program would sample every so often and compare this sample to the previous references to see which match. I don't know what the alarm recordings would be, but I'm assuming that these repeat after a period. For the comparison method, I've used cross-correlation in the frequency domain and it seems to work fine.
I've also looked at Shazam's fingerprinting algorithm (uses the peaks of the signal's spectrogram as the fingerprint) but I'm not sure if I should use it or not. I don't have much experience in signal processing and it's a complex algorithm for me, and some of the advantages it offers I don't find that important. It can compute the correlation between a large number of signals in a very short space of time (I'm not expecting a large amount of references) and it can carry out the comparison using a short amount of recording time, which I guess would be useful.
What I'd like to know is what method(s) you would recommend.
Thanks!


texts: Why are there so many windowing functions?

texts: Many windowing functions are listed here in the Mathematica documentation. I tried using a few to reduce leakage when computing a Discrete Fourier Transform. From what I could tell it made little difference which windowing function is used. Two of them that are very dissimilar are the BartlettHannWindow and BlackmanHarrisWindow. Could someone provide an example where the BarttletHannWindow is a very good choice, and another where the BlackmanHarrisWindow is a very good choice. That might help me understand why there are so many to choose from.


texts: how can I analyze a basic sine wave for its frequency phase and amplitude at a particular point in time in matlab

texts: How can I analyze a basic sine wave for its frequency phase and amplitude at a particular point in time in MATLAB? Are there any tools to do that? I am trying to do additive synthesis of inharmonic sounds in MATLAB, but having a great difficulty to extract the needed phase amplitude and frequency data of the partials I am trying to generate. Can someone shed some light regarding signal analysis in MATLAB?


texts: What is the exact compression that is happening from the first sequence of raw audio waveform data to the second sequence?

texts:
The first sequence of audio samples is recorded in JavaScript from the browser in mono, PCM format, 32-bit, and 96,000 Hz.
I wrote this audio file as a blob to a server through a JavaScript FormData object using ajax.
When I retrieved the audio from the server, I received the second sequence of audio samples. It has been downsampled to 48000 Hz and altered somehow.

Raw audio samples:


texts: Why do we assume zero mean noise in sensor data?

texts:
I am reading a paper on measuring respiratory patterns from video data. In defining the model, the authors formulate the problem mathematically as:

  $x_i(t)=h_i(t) \ast g(t) + n_i(t) $

Where $n_i(t)$ is a zero mean noise. My question is why do we model the pixel data as having zero mean noise? What prevents it from having a mean?


texts: A question about Gibbs phenomenon example in Richard Lyons book

texts:
If anyone has a copy of this book could you please shed some light here: I was trying to reproduce the example shown in Chapter 5, section 5.3.1 Window Design Method's Figure 5-19 which illustrates the ripples in the frequency domain in octave. Just as the book suggested, when you have any dropped tap in time domain, you will see Gibbs phenomenon in the frequency domain. However, my question is that, why don't you use all 32 taps in this case (the book is dropping 1 tap from 32 taps to get 5-19(c) )?  When I use all 32 taps, I can see the un-distorted frequency domain here since the iDFT and DFT pretty much cancelled out each other here...I think I am missing some points here the book is trying to make...

Thanks,

Will


texts: How do I implement cross-correlation to prove two audio files are similar?

texts:
I have to do cross correlation of two audio file to prove they are similar.
I have taken the FFT of the two audio files and have their power spectrum values in separate arrays.

How should I proceed further to cross-correlate them and prove that they're similar? Is there a better way to do it? Any basic ideas will be helpful for me to learn and apply it.


texts: Differential BPSK vs BPSK what are the advantages and disadvantages?

texts:
Recently I have come across different variation of the sane modulation techniques. I am quite confused about the advantages of different constellations. In particular, what are the advantages of using

Differential BPSK over BPSK?

Any thoughts?


texts: Equation for impulse train as sum of complex exponentials

texts:
Could someone please break down what's going on in this equation for me? I understand what the left side looks like, but not so much how the right side is the same thing.

Impulse train:

$$\sum_{m=-\infty}^\infty \delta[n-Nm]=\frac 1N\sum_{k=k_0}^{k_0+N-1}e^{ \ j 2\pi kn/N} \\ \\ n, m, N,k_0 \in \mathbb{Z}$$

where $\delta[n]$ is the Kronecker delta:

$$ \delta[n] = \begin{cases}
    1 \quad n = 0 \\
    0 \quad n \ne 0 \\
\end{cases} \\
\\
n \in \mathbb{Z}$$


texts: Don't understand the relation about the wireless communication flow chart and TDMA illustration diagram

texts: The first picture is the wireless communication flow chart,and the second is the TDMA illustration diagram,Which part of first picture is TDMA ? or if first picture is the standard process for TDMA,then how to explain it by using the TDMA illustration diagram? I mean, TDMA illustration diagram just show me  that divide whole time into several part ,and it seems that there is no relation about the block in the first picture like source coding,modulation ,etc.


texts:
Say a narrow band signal $n(t)$ has Power Spectral Density (PSD) $S(f)$.

If the signal $n(t)$ got multiplied by $\cos(2\pi Ft)$, then what will the PSD of resulting signal in terms of $S(f)$ be? Will this be $\frac{S(f+F) + S(f-F)}{2}$  ?


texts: Plotting a Time Frequency contour/colormesh plot of a Discrete Wavelet Transform

texts:
I have a pressure vs time data of a noise on which I wish to perform discrete wavelet transform.
I have divided my frequency range into 1/3rd Octave Bands and have calculated sound pressure level at each band.
I am very much confused on how to perform a discrete wavelet transform in 1/3rd Octave Frequency Bands or any set of bands for that matter and further how to interpret the output(coefficients) so as to plot a Time-Frequency colormesh or a contour plot of my results.
I was able to plot a continuous wavelet transform of my data using the pcolormesh module of matplotlib. The output of pywt.cwt are easy to interpret and hence I could plot a colormesh for it.
I am a beginner in DSP and have just started learning wavelet transforms, so any help will be much appreciated. Thankyou.


texts: YIN pitch estimation algoritm-simplified explanation

texts:
can someone please explain the steps of YIN pitch detection algorithm in a simple way "especially the last 3 steps".

here is the research paper of YIN algorithm
http://audition.ens.fr/adc/pdf/2002_JASA_YIN.pdf
Thanks in advance

Edit1(as i understand):
I understand the autocorrelation to find out if the signal periodic or not "step1", then we will take the difference between the original signal with zero lag"delay" and subtract it with a delayed one "step2",now we have the negative peaks but to make the algorithm works and pick up the first delayed negative peak with the lowest delay we have to remove the first negative peak with 0 delay"step3", the fourth step is to set a threshold to get the most negative peak(i am not sure if i am right)"step4",step 5 & step 6 i dont understand them,can you please tell me if i am explaining right or wrong


texts:
I have the signal $s(t) = \cos(2\pi f_1t)$ and I am looking for its components vs the Fourier basis, over the interval $[0, T]$. The formula for computing the coefficients is
$$
s_n = \int_{t_0}^{t_1} s(t) \frac{e^{-j\frac{2\pi}{T} nt}}{\sqrt{T}} dt
$$
and in this case $t_0 = 0$ and $t_1 = T$. The solution seems to be
$$
\frac{\sqrt{T}}{2}(-1)^n\bigg[\mathrm{sinc}\big(f_1T+n\big)e^{-j\pi f_1T} + \mathrm{sinc}\big(f_1T-n\big)e^{j\pi f_1T}\bigg]
$$
but i can get no further than this
$$\frac{\sqrt{T}}{2}\bigg(\frac{e^{j2\pi (f_1T-n)}-1}{j2\pi (Tf_1-n)} + \frac{e^{-j2\pi (f_1T+n)}-1}{-j2\pi (Tf_1+n)}\bigg)$$


texts:
Assume, I am trying to use a 50khz sampling rate for some signal. I don't know what the spectrum of input signal is, so I will have to low pass manually with sinc (theoretical best low pass) to band limit it to 25khz.
Ignoring noise for now, assume my input signal is a 5khz signal starting at time = 5s, all values before that 5s is 0. Now to sample this, I would theoretically be passing this over a sinc low pass filter at 25khz and then sampling right?
At steady state (time = infinity) when I apply this low pass my output will just be the same sine at 5khz but with same amplitude. At steady state, no matter how much time I apply this low pass I'll still get a 5khz sine with same amplitude, because I have a train of 5khz sines preceding and succeeding it.
But during transients, ie at 5s areas I will be getting this sine smeared due to low pass right? In these transient scenarios, I am thinking every instance of such low passing will smear it further.
Would convolving a sinc used for band limiting at 25khz with itself would give me a result that explains this smearing? An ideal infinite bandwidth sinc would just have it's time domain property as a rectangular pulse, so multiplying them with each other will still return a rectangular pulse. Convolving two identical infinite bandwidth sincs will give back a result that's the same as the original sinc. Bandwidth is infinite and all transients are preserved, so it makes sense.  In fact I'm getting back just the original signal.
But for a band limited sinc, the time domain property would also have Gibbs phenomenon on top of the rectangular pulse right, and if I multiply it with itself I'll actually get it even more smeared is my understanding. If I convolve two identical band limited sinc (limited at say 25khz) would the result be the same as original or would it be more smeared?
Practical scenarios of non infinite time is likely to have even more deviations, but for now I would like to know how the response would be for the kind of signal and process described above.
Note: there is no jump discontinuity. The sine just starts from time 5s, the value at 5s is 0. Till 5s the whole signal is a dc value of 0, and after 5 second it becomes a 5khz sine having value 0 at 5seconds.


texts: Inverse fourier transform by using a random phase and an original magnitude (matlab)

texts:
I want a image. Firstly I want to do fourier transform and get the magnitude and phase information. And then make a new random phase, finally use the changed random phase and original magnitude to make the inverse fourier transform. But I can not achieve it?

The code is shown here.

 img   = imread('flag.jpg');
 imafft = fftn(double(img));
 mag = abs(imafft);
 phi= angle(imafft);
 s= size (img);
 rand_phase= 4*pi*(rand(s(1),s(2))-0.5);
 newfft = mag.*exp(i*rand_phase);
 ima_out = ifft2(newfft);
 figure(2);
 imagesc(abs(ima_out));
 axis image;
 title('Ramdom phase version');

The problem is in the sentence of newfft = mag.*exp(i*rand_phase); I do not know how to modify.

Who can help me? Thank you so much.


texts: Stable digital filter with leading zeroes

texts:
I'm looking construct a stable pole-only filter where the feedback coefficients start with a block of zeroes, i.e.

\begin{align}
a_0 &= 1\\
a_i &= 0, \textrm{ for}\ 1 \le i \lt k\\
a_i &\in\mathbb{R}, \textrm{ for}\ k \le i \le n
\end{align}

For $k = n$, this is a single-tap feedback delay, e.g.

$$
y_i = x_i + a_k y_{i - k}
$$

which is stable as long as $|a_k| < 1$.

However, I'm looking for a way to construct such filters with arbitrary $k$ and $n$.  Are there conditions or constraints that would let me construct this?


texts: How can I determine the fundamental period of the signals?

texts:
How can I find the fundamental period of the given signal? Can anyone help me, I am stuck.

$x(t) = cos(2πt) + cos(4πt)$
$y(t) = cos(2πt) + cos(2t)$

Thank you in advance.


texts: Which modulation technique is better for long distance transmission, if I'm using the CC1101?

texts:
I'm using the Texas Instrument transceiver cc1101 at 430-440 [MHz] and I'm wondering what the best modulation technique is for long distance transmission.
This transceiver offers several modulation schemes:
2-FSK
4-FSK
GFSK
MSK
ASK
Ook
What should I pay attention to if I want to use the best one to transmit for at least 5 [Km]?
Extra information:
I am using a power amplifier to achieve increased transmission distance, I also know that a lot falls into the aspect of using directional antennas and other factors such as direct line of sight. Knowing this, I want to know if what I should focus on to achieve my goal is the BER of these techniques or something else.


texts:
In Alaa Kharbouch, Ali Shoeb, John Guttag, Sydney S. Cash,
An algorithm for seizure onset detection using intracranial EEG,
Epilepsy & Behavior,
Volume 22, Supplement 1,
2011 (section 2.1, 3rd paragraph) about EEG, the authors note that the spectral amplitude profile of a signal is inversely proportional to frequency. To correct for this trend, they propose to apply a derivative filter to the signal.

My question is: Is there a Python function which implements such a derivative filter ? Is the savgol_filter function from the Scipy module suited to this task ? If not, how could I design such a filter in Python ?


texts: Analog butterworth to digital - bilinear transform

texts:
In my previous question I've designed analog butterworth filter (poles own calculated). But now I would like to transform it to digital domain. I'm using bilinear transformation but not all is clear for me. Please look at my script below:

clc, clear all, close all;

%ANALOG BUTTERWORTH FILTER DESING
%https://engineering.purdue.edu/~ee538/AnalogFilterDesigns.pdf

%Butterworth LP with -3dB at fc and N order and G gain

wc = 1000;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain

%analog filter design with prewarped frequency
%poles and zeros

%frequency normalizing and prewarping
wc = 2*tan(wc/2);

%poles searching
for k = 0:N-1
   sk(k+1) = wc*exp(j*(pi/2))*exp(j*(2*k+1)*(pi/(2*N)));
end

[NUM,DEN] = zp2tf([],sk,G);
b0 = DEN(end)*G      %gain normalization

Hs = tf(b0,DEN)

figure();
freqs(b0,DEN)
title('lowpass filter with pre-warped frequency visualization')
figure();
plot(sk./10^2, 'x')
xlim([-10,10]);
title('Poles placement in s-plane')

%continuous to digital conversion
[NUMD,DEND] = bilinear(b0,DEN,1);

fvtool(NUMD,DEND)

Results looks like below:

I'm sure that is not what I want, because I have nowhere applied frequency sampling that is interesting me. I'm not also sure if I've done prewarping well.

EDIT: frequency response for ws = 2000 in bilinear transformation:


texts: If $\theta$ is uniformly distributed in $(0, 2\pi),$ then what is the distribution of $e^{i\theta},$ where $i = \sqrt{-1}?$ And what are the statistical properties of $\left[e^{i0\theta}\, e^{i1\theta}\, e^{i2\theta}\dots\, e^{i(N-1)\theta}\right],$ where $\theta$ is a single observation and $N$ some integer?


texts: Deconvolution of non-stationary, 1-D signal?

texts: I have a time series that has been measured after convolution with a moving average filter. Knowing the parameters of the moving average filter, is it possible to reconstruct/constrain the values of the original? When considering deconvolution, there is much discussion in the literature of the form and magnitude of measurement error, but the variation in signal due to this component is much less than the non-stationary component (generated by the physical process generating the signal), if that knowledge helps in the analysis. This is not my primary domain of expertise so my apologies if it is a simple question (even getting a gauge of whether this is simple or difficult would be helpful). I wonder if I need to look into Gabor deconvolution or other topics that would help me understand an approach to this problem. Thanks very much.


texts: Determining Fourier Series coefficient for Discrete time

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

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

Please, anyone, help me to understand this


texts: how does the Frequency spectrum of a signal change in the logarithmic space?

texts:
I am trying to implement basic homomorphic filtering but I can't seem to understand what happens in the frequency domain when you take the natural logarithm of the function.

In Matlab I'm getting not a number or infinity everywhere.


texts: How to uniquely identify phase difference between harmonics?

texts:
Let's say I have a multi-frequencies signal of 440Hz and 880Hz, slightly out of phase.

I can take a chunk of the signal and calculate phase difference of the two frequencies with FFT bins. However, since frequencies are different, the phase difference will differ depending upon where I took the chunks, despite they are taken from the same signal.

I wonder whether there is a way to uniquely identify phase difference between harmonics, with arbitrary chunks? (so that I can identify harmonic signals.)


texts: Simulate BER curves for OFDM with convolutional coding

texts:
I am tasked with generating performance curves for an OFDM transmitter/receiver, and I am unsure how to calculate the appropriate AWGN power, and where exactly the apply it.

The transmitted packets contains typical non-payload portions, such as a synchronization marker, short training sequence, long training sequence, cyclic prefix, pilot subcarriers, and guard band subcarriers. The models have quite a bit happening in terms of real-world effects: Carrier Frequency Offset simulation, multipath simulation, channel-specific effects simulation.  These are used to exercise the receivers ability to handle these real world effects, using CFO correction as well as Equalization.

Furthermore, the payload bits are applied to a 1/2 convolutional encoder/decoder, so the IQ symbols look like QAM-4 (2-bits per subcarrier), but after the decoding happens, it's an effective 1-bit per subcarrier.

How can I apply noise the transmitted waveform to correctly generate a BER vs EbNo plot?


texts: Calculating the first derivative of an image using DFT

texts:
I need to calculate the first derivative of a greyscale image (a 2D array) using a DFT function I built (which works). Unfortunately, the results don't seem to be correct.

In the fourier domain, the derivative d/dx is given as F(u,v)* 2*pi*i/N * u, where u is the x-axis transformed, N is the size of one of the matrix's dimensions, v being the other one.

Attached is the code. What bothers me is that I'm not getting the same results as I would differentiating by convolution with (1,-1) or (1,-1) as a column vector.

def derivative(fourier_signal):
"""
Derivative in fourier domain is multiplying by u or v, and 2pi*i/N
:param fourier_signal:
:return:
"""
N = np.shape(fourier_signal)[ZERO]
M = np.shape(fourier_signal)[ONE]
u = np.arange(N)
v = np.arange(M)
du = fourier_signal * (u*TWO_PI*1j)/N
dv = fourier_signal * (v*TWO_PI*1j)/M
return du, dv

def fourier_der(im):
    # Calculate DFT2
    dft_image = DFT2(im)
    # Function that Multiply by rows by u, columns by y
    du, dv = derivative(dft_image)
    shifted_du, shifted_dv = np.fft.fftshift(du), np.fft.fftshift(dv)
    dx, dy = IDFT2(shifted_du), IDFT2(shifted_dv)

I'm not looking for easy answers on how to do it, but rather a direction to as to why my output is incorrect.


texts:
I'm using Matlab to process neurophysiological data and i'm currently dealing with this sort of noise in my signal:

Close up

EDIT:

For clarification, the clean signal should look like this:

EDIT 2:

Clean

Noise contaminated

These spectra were taken from 10s segments as the whole block is considerably bigger. Both spectra are from segments that include the data in the first and third figures (with and without noise bits of data). The noise events are sparse within each data block.

This happens somewhat frequently across datasets and due to it, i'm having some trouble detecting events of interest (i.e. neural spike detection). As i'm trying to detect small amplitude spikes in the signal, it's critical to preserve as much of the original signal as possible.

Any suggestion on how to better approach this problem? What would be te best way to correct these shifts/spikes in the signal, keeping as much of the signal as possible?

Any suggestion/help will be immensely appreciated!

Cheers, M


texts:
What is the difference between implementation a Reed-Solomon code/encoder/decoder accorsing CCSDS standard and DVB?
I am not sure that the differences are the lenght of a codewords and rate.


texts:
I have a signal with length of $N$.
The sampling frequency is $f_s$.
Window size is $w$
The noverlap size is set to $m$.
And $nfft$ is selected for frequency resolutions.

I am running a STFT with matlab spectrogram function.

[S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)

what is the effect of nfft , $f_s$ and  noverlap on STFT
computational burden?
what is computation complexity of STFT with these parameters size?
{in term of big O() }

Note: The FFT computational complexity for a signal with length of $N$ is $Nlog(N)$


texts:
I have a few questions about this assignment I am solving.
It is q3.33(a) from oppenheim's DSP book 2nd edition.

basically, the question gives you a pole zero plot, from which I was able to determine the correct x[n]. (Couple classmates confirmed my answer on that).

and it also gives you the y[n] of the form $(0.5)^n \cdot x[n]$.

And wants you to work out a pole-zero of y[n] in the z plane.

my question is this: given that I already know X[z], can i just do z transform of $0.5^n$ as if it was $a^n \cdot u[n]$?

Then I will have $[\frac{1}{1-0.5z^{-1}}] \cdot X[z]$

Thanks everyone.


texts:
I'm going over the W3 PNG specification (creating a PNG library from scratch) and I've finally found how to create a green 1x1 image.

Now I'm trying to create a bigger image of mixed Red, Green and Blue pixels. Let's say a 4x4 image. Sadly I'm getting all the pixels mixed and some of them are Black or Pink.

Details:

Signature: OK
IHDR: OK

Width: 4
Heigh: 4
Bit depth: 8
Color: 2
Filter: 0
Compression: 0
Interlace: 0

IDATA: 4 chunks*
IEND: OK

*IDATA chunks:

4 bytes length: Number of bytes of the DATA after Zlib deflate
4*1 byte type: IDAT
X bytes data: 4*4bytes unsigned integers, Zlib deflated, one after the other, Network byte order

255 for blue**
65280 for green**
16711680 for red**

4 bytes crc: OK

**Bitwise results of:

alpha<<24 | red<<16 | green<<8 | blue

With alpha, blue, green and red taking values from 0 to 255

What's wrong with this?


texts: Redundancy when using Nyquist rate for a time series

texts: I sample temperature at one sample per second (the hardware I am using takes the temperature at this rate, so this is the max I can efficiently sample.)  The Nyquist rate for this signal would be 2 samples / second, correct?  If so, both samples (for any given second) will be the same, and so redundant.  For discrete time series does the Nyquist rate apply?  I am sure my logic is wrong somewhere, but what am I missing?


texts: How to design a continuously time-variant digital filter?

texts:
I have discrete time series containing signal with smoothly varied frequency over time (called a "sweep"). How can I design a discrete filter (low-pass or band-pass in my case) of a finite length with linearly varying cut-frequency over time and constant cut-slope?

EDIT: the signal is the sampled "trace" of the seismic source - a seismic vibrator, which sends the vibrations of the slowly varying frequency down the earth. The dependency of the frequency over time (the sweep) is known (let it be linear, $f(t)=f_1*(1-t)+f_2*t)$, but there is a problem that there might be another vibrators that operate on their own, and the task is to "band-guard" the trace of this vibrator avoiding the unwanted signals from other ones.


texts: Why does it make sense to have noise of a different shape than the state vector?

texts:
Usually when I see the Kalman filter in lectures (or other resources like 1), the state equation is

$$\mathbf{x}_{k+1} = A_k \mathbf{x}_k + B_k a_k + r_k^{(s)},$$

and the measurement equation is

$$z_k = H \cdot \mathbf{x}_k + r_k^{(m)}$$

where $r_k^{(m)}, r_k^{(s)}$ is normal distributed noise of the same shape as the state vector $\mathbb{x}_k$.

Then, the covariance in the prediction step is usually stated as

$$P_{k+1}^{(P)} = A P_k A^T + C_k^{(r_s)}$$

where $C_k^{(r_s)}$ is the covariance matrix of the system noise.

However, in one lecture the state equation was stated as

$$\mathbf{x}_{k+1} = A_k \mathbf{x}_k + B_k a_k + G r_k^{(s)}$$
and the prediction for the covariance matrix of the system noise was
$$P_{k+1}^{(P)} = A P_k A^T + G Q G^T$$

Why does it make sense to model the noise as something which gets multiplied by a matrix $G$ in contrast to simply having a noise vector which has the same shape as the state vector?

(See my article for detailed explanations of the variables / how I understand the Kalman filter right now.)


texts: Is it normal to see periodic peaks going from 48.75 - 90 MHz? Also i was checking on my SDR dongle once and the whole band was raised above background offset for like a whole 50 MHz bandwidth and then it dropped to background again? Does anybody know what that means?


texts:
I'm successfully soft-decoding D-BPSK by taking the dot-product of the constellation-position of the symbol and of the previous symbol.  If the result is >= 1, then the symbol phase hasn't changed and the bit is a zero.  If the result is <= -1 then the phase has shifted and the result is a one.  In between -1 and 1 the result is a soft 0 or soft 1.

I can't figure out how to do the same thing with D-QPSK.  I can use just the phase, but this throws away a lot of information that could help the soft-decoder.

This paper explains how to do it and gives a formula (10):

$b_1 = \mathrm{Re}\{s_n s^*_{n-1}\}, b_2 = \mathrm{Im}\{s_n s^*_{n-1}\}$

But I don't understand the notation — what does a * floating above mean?  I tried just multiplying the complex numbers and taking the real and imaginary parts but this didn't work.

Since the constellation can rotate, how can the two axes be teased apart?


texts: Non-orthogonal basis which offers better energy compaction in image compression

texts:
I am aware of image-independent basis, i.e. DCT, and image-dependent basis, i.e. Karhunen–Loève, which are used in compacting energies for image compression. These bases are orthorgonal.

Are there any compression basis which are non-orthogonal that delivers better compaction than the optimal KL transformation? They can either be independent or dependent upon the statistics of the image.

EDIT: I have forgotten to mention that I am looking at block-based transformation coding.

The motivation of the question is due to the following paper Optimized Nonorthogonal Transforms for Image Compression. The premise of this paper is that the block based Karhunen–Loève transformation yields compacted coefficients, but the coefficients are still correlated across different blocks as it is often the case where adjacent blocks of a natural image are correlated to one another. To eliminate correlation across blocks, the optimal transformation would involve the use of non-orthogonal basis.


texts: Generating sound of particular frequency matlab

texts:
I want to generate a high frequency (single) tone. I am using matlab (any other platform/tool for the same will also work). I am using the following snippet which I wrote by referring to various sources achieving similar specification.

amp = 10;
Fs = 44100;
duration = 60;
freq = 20000;
t = (0:1/Fs:duration-1/Fs);
%t = linspace(0,duration,Fs*duration);

wave = amp*sin(2*pi*freq*t);
%sound(wave, Fs);

audiowrite('20KHz-pilot.wav',wave,Fs);

draw_spec = 0;
if(draw_spec == 1)
    NFFT = 2048;
    wave = wave(1:NFFT);
    y = fft(wave);
    y = abs(y/NFFT);
    y = y(1:NFFT/2+1);

    % as 0 and Fs don't have any alias while others do
    % we need to double the power of the other frequencies since we are
    % converting from 2 sided power spectrum to single sided spectrum.
    y(2:end-1) = 2*y(2:end-1);

    f = Fs*(0:(NFFT/2))/NFFT;
    plot(f,y);
end

The problem I am facing is that when I load the file in audacity it shows frequencies other than 20KHz in fact 20Khz is not there in spectrogram. Further when I play the sound I can hear it which should not be the case (most adults have audio hearing range upto 18Khz). Also plotting spectrum of the wave in matlab only shows 20KHz.

I have tried to use common conventional variable names so the code should be pretty self explanatory.

I am not able to find where I am making a mistake, is there a mistake in the code or is it that due to sampling rate etc. there are other frequencies also generated (aliasing or similar concepts). I tried playing the sound on different speaker systems to check if it could be a hardware issue but that was not the case.

NOTE: I have browsed through multiple resources with similar question title as mine, unfortunately none of them were able to help me. Thus, I resorted to posting this question.


texts: DFT of real sinusoids - why sum over -$N/2$ to $N/2-1$ as opposed to $0$ to $N-1$?

texts:
I'm going through a Coursera course on signal processing, and we're just introduced to DFTs.

We are told that if you have a complex sinusoidal signal
$x[n]$ where $n=0,1...\ N-1$, its DFT is given as

$$X[k]=\sum_{n=0}^{N-1} x[n] e^{-j2\pi kn/N}$$

And this makes sense to me. Because $x[n]$ has $N$ samples. You take each one, multiply it with the sample of a basis function, essentially taking the dot product of the signal with the basis function, and you get how much of that frequency is present in the input signal.

But we are told that if you have a real sinusoidal signal $y[n]$ where $n=0,1...\ N-1$, then the DFT is given as

$$X[k]=\sum_{n=-N/2}^{N/2-1} y[n] e^{-j2\pi kn/N}$$

My question: Why do we sum over the different values? What is the intuition and explanation behind this?

I'm not even sure what $y[-N/2]$ would mean.

Is the $-N/2^{th}$ sample the same as the $N/2^{th}$ sample because we assume the input signal is periodic? Even if so, why is taking the DFT summation over such a range beneficial? We could also do this for complex sinusoidal signals since they're periodic too right? What's the advantage?


texts:
Prove that the convolution of two even functions is an even function.

I have my own proof which I have included as an answer, but it assumes a linear time-invariant system. I want to know if there is a general way to prove this (or if it even holds for a general convolution) ?


texts: fft - spectral leakage from noise into target frequency?

texts:
I have an acoustic setup. A sender, emitting a clean sine signal. The receiver is far away. In the time domain, the sine is hardly visible because the noise is nearly as big as the signal. In the frequency domain, the noise is distributed across all frequencies.

I need to extract the amplitude at the target frequency, I do not care about the phase.
Currently, I calculate the length of my fft-window to contain exactly N periods of my target frequency (n_bins=N * samplerate / frequency) to avoid a loss in magnitude due to spectral leakage.
To remove the noise at this frequency, I do the same with a measurement w/o signal to estimate the noise background and subtract it.

My question now is, should I still use a window-function like blackman? From my understanding, I could now have leakage from noise in adjacent bins into my signal?

I have tried using a window-function in a python-example for this, with simulated data. Here, when using a window, the amplitude at my target frequency was actually higher than expected, while the results without the window were closer to reality.

EDIT:
Python source code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import blackman as w_function

target_freq = 12000
samplerate = 1000000.
window = True

# generate data
x = np.array(range(0,17166)) / samplerate
y = 290*np.sin(2*np.pi*target_freq*x) # signal
y += 150*np.sin(2*np.pi*(target_freq - 451)*x) # noise near the signal
y += 232*np.sin(2*np.pi*200185*x) # noise far away

# find closest multiple
reduced_len = len(y)
while not float(reduced_len * target_freq / samplerate).is_integer():
    reduced_len -= 1

def calc_fft(window):
    # prepare window
    if window:
        wndw = w_function(len(y))
    else:
        wndw = np.array([1] * len(y))
    wndw_correction = float(len(y)) / float(np.sum(wndw))  # amplitude     correction factor for window

    # calculate fft, apply corrections for correct amplitude     representation
    rft = np.abs(np.fft.fft(y * wndw, reduced_len, norm=None))
    rft *= 2 * wndw_correction / float(len(rft))

    # calculate corresponding frequencies, show only positive frequencies
    fft_freq = np.fft.fftfreq(len(rft), 1 / samplerate)
    rft = rft[np.argsort(fft_freq)]
    fft_freq = fft_freq[np.argsort(fft_freq)]
    low_ind = (np.abs(fft_freq-0)).argmin()
    fft_freq = fft_freq[low_ind:]
    rft = rft[low_ind:]
    return fft_freq, rft

fft_freq, rft_window = calc_fft(True)
fft_freq, rft_no_window = calc_fft(False)

plt.plot(fft_freq, rft_window, label="window")
plt.plot(fft_freq, rft_no_window, "-.", label="no window")
plt.legend()
plt.show()

Without any noise (just comment lines in python out), the fft without window gives the exact amplitude (290), with window the amplitude is to high:

With noise, the amplitude with window is still to high, and the amplitude without window is somewhat to low:


texts: How to choose a phase for the deconvolution of an autocorrelation?

texts:
Say I have a function, $C=C\left(x\right)$, whose fourier transform is denoted by $c=c\left(k\right)$, i.e. $C\left(x\right)=\sum_{k=-\infty}^{\infty}c\left(k\right)\chi\left(x\right)$, where $\chi\left(x\right)$ is some complex valued basis function.

I know that the autocorrelation is given by $M\left(\Delta x\right)=\sum_{k=-\infty}^{\infty}m\left(k\right)\chi\left(\Delta x\right)$, where $\Delta x = x_2-x_1$ and by the convolution theorem I know that $m\left(k\right)=\left(c\left(k\right)\right)^*c\left(k\right)$. This means that knowing the coefficients of the autocorrelation immediately gives me the squared magnitude of the original coefficients. This implies that any two functions $C'\left(x\right)=C\left(x-x'\right)$, that differ only by a translation, will possess the same autocorrelation, $M\left(\Delta x\right)$, and so if I am trying to obtain the $c\left(k\right)$ from the $m\left(k\right)$, then I can choose an arbitrary phase for the $c\left(k\right)$. However, I don't think I can choose the phase of $c\left(k\right)$ independently for each $k$, it seems like there needs to be some way to consistently choose the phase for all $k$ (which corresponds to specifying a translation for $C\left(x\right)$).

My question is: How do I choose a consistent phase for each of the coefficients for the original function, i.e. for all $k$?


texts: Resources on Solving Convex Optimization Problems in the Compress Sensing Field

texts:
When I read papers of compressed sensing, sparse representation and whatever requiring optimization of a cost function, I just find the final results as an iterative equation or so which will converge after few iterations and solve the problem.

But, when I try to add some other constraints to the cost function, I cannot derive the iterative or whatever the solution is and I know that is for my poor mathematics.

Can anyone lead me to a book or so to read and understand how to solve these types of mathematics? What book chapter? What tutorial?


texts: NCO loses lock after a while for certain values of alpha

texts:
EDIT: So the transfer function for my system is

$$ F_{nco}(z)=\beta E(z)/(z-1)$$
$$ z\Theta_{nco} (z)=\Theta_{nco} (z)+\alpha E(z)+\beta E(z)/(z-1)$$

So I would like to obtain $\Theta_{nco} (z) / \Theta(z)$ transfer function cause I think that would allow me to calculate whic range of $\alpha$ would be good for my PLL to lock. The thing is I don't really know how to obtain $ E(z)$ from $e[n]= \arg \big\{ x[n] \cdot (y[n])^* \big\} $

I'm trying to code my own NCO in order to use it for my own Costas loop (since I couldn't get this to work Implementing Costas loop for 4-PAM) and so far the NCO is working good except for it losing lock after a while. These are the update equations I'm using:

$$\theta_{nco}[n+1]= \theta_{nco}[n] + \alpha e[n] + f_{nco}[n]$$

$$f_{nco}[n+1] = f_{nco}[n] + \beta e[n]$$

So I generated this input with $f=50 Hz$ and $f_s = 1 kHz$ during $0.2 s$ (that is, 200 samples):

$$\begin{align}
x[n] &= e^{j 2\pi n f/f_s+\theta} \\
  &= \cos(2\pi n f/f_s+\theta)+j\sin(2\pi n f/f_s+\theta) \\
\end{align}$$

Then I'm running a for loop to update the NCO, which initially has 0 phase and frequency ($\theta_{nco}$ and $f_{nco}$). The parameters being $\alpha = 0.1$ and $\beta = 0.0025$.

The way I obtain the error is

$$e[n]= \arg \big\{ x[n] \cdot (y[n])^* \big\} $$

being $y[n]$ the output of the NCO:

$$y[n]=e^{j (2 \pi n f_{nco}[n] + \theta_{nco}[n])}$$

The thing is, it locks quite fast but it loses lock towards the few last samples (at around 0.2 s) (see attached picture). If I increase the number of samples, it loses lock even more.

I managed to fix it by changing alpha and beta to $\alpha = 0.05$ and $\beta = 0.000625$

So basically what I'd like to know is how I can know which are the best values for $\alpha$ and $\beta$ since I'm just randomly picking a low number for $\alpha$ and then doing $\beta=\alpha ^2 /4$. Thanks!


texts:
I currently have a real-time noise suppression algorithm based on two inputs from different microphones. It has passed most tests that I have thrown at it, however when the signal I am trying to retrieve has a slight delay one mic relative to the other it works poorly.

I have done preprocessing delay estimation in the input signal through cross-correlation just for testing and it works, however in the final algorithm I can't preprocess the signal as it is supposed to be real-time.

I am taking sub windows of signals and applying my algorithm, what I tried to do is apply the cross-correlation to each sub window, it takes the FFT of the sub windows does the calculation and then the IFFT just to get the delay each time. I later have to take the FFT of the sub window with a different windowing function for my filtering algorithm. This much extra processing will drastically reduce the sub window size my algorithm can keep up with in real-time.

The question is, is there a faster way to do this delay estimation? Maybe something that works with each input sample to get a better estimation each time? Something like a control loop?

Thanks in advance!


texts: relationship between number of wavelet decomposition level and noise level

texts:
I am studying wavelet theory by myself and hope to use wavelet to denoise images. I write one matlab program to watch denoised image.

Gray = ones(512, 512, 'uint8');
Gray = Gray.*128;
Gray = imnoise(Gray, 'gaussian', 0, 0.01);
[cA1, cH1, cV1, cD1] = dwt2(Gray,  'db1');
[cA2, cH2, cV2, cD2] = dwt2(cA1/2, 'db1');
[cA3, cH3, cV3, cD3] = dwt2(cA2/2, 'db1');
[cA4, cH4, cV4, cD4] = dwt2(cA3/2, 'db1');
[cA5, cH5, cV5, cD5] = dwt2(cA4/2, 'db1');
figure; imshow(cA1/2, [0, 255]);
figure; ismhow(cA2/2, [0, 255]);
figure; ismhow(cA3/2, [0, 255]);
figure; ismhow(cA4/2, [0, 255]);
figure; ismhow(cA5/2, [0, 255]);

I find cA3 contains a lot of noise even if three decompositions have been performed. If I increase noise level, I need to perform more decomposition to obtain one smooth image in the approximation coefficients. If number of wavelet decomposition level is big, modifications of detail coefficients of high decomposition level will affect a lot of pixels and edges may be smoothed. I try to change wavelet family,but the problem still occurs.
By the way, the function dwt2 seems to magnify the values of approximation coefficients by 2. Why does dwt2 behave like that?


texts: Why does BLE require twice the bandwidth of Bluetooth Classic?

texts:
Bluetooth Classic Basic Rate sends data at a rate of 1 Mbps using GFSK modulation in 1 MHz-wide data channels. Bluetooth Low Energy also transmits at that same data rate and modulation, however, it requires 2 MHz-wide channels. Why?

Is this about providing less stringent specifications so cheaper filters can be used in the design?
Is it about providing a better SNR?
Lastly, does the information actually use the full bandwidth or could I still see the data through a 1 MHz filter centered on the channel albeit at a reduced SNR?


texts: derivation of discrete Fourier Transform

texts:
I'm learning DSP on my own. Found a free online book. Currently reading a chapter on Discrete Fourier Transform. Right out of the gate I can't understand several things. I'm pretty good with math and have a decent understanding of complex numbers and their representation with the help of $e^{j\phi}$. But this got me stumbled.

Below is the paragraph I can't understand and the questions pertaining to it.

  We now develop a Fourier representation for finite-length signals; to
  do so, we need to find a set of oscillatory signals of length $N$
  which contain a whole number of periods over their support. We start
  by considering a family of finite-length sinusoidal signals (indexed
  by an integer $k$) of the form

  $$ w_k[n]=e^{j \omega_k n}, \space n = 0,...,N-1 $$ where all the $\omega_k$’s are
  distinct frequencies which fulfill our requirements. To determine
  these frequency values, note that, in order for $w_k[n]$ to contain a
  whole number of periods over $N$ samples, it must conform to

  $$ w_k[N]=w_k[0]=1 $$ which translates to $$ (e^{j \omega_k})^N = 1 $$

  The above equation has $N$ distinct solutions which are the $N$ roots
  of unity $e^{j\frac{2\pi m}{N}}$, $m = 0,…,N - 1$

which contain a whole number of periods over their support.
What does "periods of their support" mean?
note that, in order for $w_k[n]$ to contain a whole number of periods over $N$ samples. Why does $w_k[n]$ must contain a whole number of periods over $N$ samples? Where does this requirement come from?
it must conform to $w_k[N]=w_k[0]=1$. Why? Where does this requirement come from? I do understand $w_k[N]=w_k[0]$, that's because it is periodic, what I don't understand is $=1$ portion. Why at sample $0$ and sample $N$ it must be equal to 1?
which translates to $(e^{j \omega_k})^N = 1$. How does $w_k[N]=w_k[0]=1$ translate to $(e^{j \omega_k})^N = 1$?
The above equation has $N$ distinct solutions which are the $N$ roots of unity $e^{j\frac{2\pi m}{N}}$, $m = 0,…,N - 1$. Can someone demonstrate the finding of these roots?


texts: A general conceptual understanding of moving average

texts:
I have a conceptual question here. I understood how to calculate a Moving Average(MA).

For example :

x(n) = 0, 1, 2, 3, 4, 5

MA = 0, 1/2, 3/2, 5/2,7/2, 9/2, where for calculation of the MA we only take the past values without any update

Is there also a case , where we use the newly calculated values and use that to update the MA for the next time instant?

x(n) = 0, 1, 2, 3, 4, 5

MA (n=0) = 0

MA (n=1) = (0 +1) /2 = 1/2

MA(n = 2) = (2 + 1/2) /2 = 5/4

MA(n = 3) = (3 + 5/4) = 17/4 and so on ?

Does this find an application anywhere?

Please help me with this understanding.


texts: I read several papers on physical layer security that explain several methods to maximize the secrecy rate. However, I didn't find a very clear definition about the secrecy rate. What is it? How is it computed?


texts:
Suppose, I have 3 kernels:

$$ \left[
\begin{array}{cc}
  a&b&c\\
  d&e&f\\
  g&h&i
\end{array}
\right] $$
$$ \left[
\begin{array}{cc}
  p&q&r\\
  s&t&u\\
  v&w&x
\end{array}
\right] $$
$$ \left[
\begin{array}{cc}
  \alpha&\beta&\gamma\\
  \delta&\epsilon&\zeta\\
  \eta&\theta&\iota
\end{array}
\right] $$

How can I create a filter bank from them?

Should I AND, or OR, or add them together?

Or, should I just apply each of them one by one to my test-image through three separate convolution operations?


texts: What are the advantages of Laplace Transform vs Fourier Transform in signal theory?


texts: Difference between Gabor filters, Gabor wavelets and Log-Gaber wavelets

texts:
I am a beginner in the field of image processing, and I am trying to understand Gabor transform but the more I read, the more I get confused.

So I wish to find someone who can help me answering this question: what is the difference between a Gabor filter bank, a Gabor wavelet and a Log-Gabor wavelet?


texts: When can the impulse response become zero?

texts:
The article Efficient Use Of Sparse Adaptive Filters (Proc. Asilomar Conference, Khong et al., 2006) introduces adaptive filters for the estimation of channels or systems having a sparse impulse response. Assuming that the channel is modeled by an FIR model, then the output of the channel will be given
by $$y[n] = \sum_{i=0}^{p-1} h[i] x[n-i] \tag{1}$$ where $x$ is the input to the model and $h$ is the MA coefficients. As an example, consider an MA(3) model of order p=3 having coefficients h = [1,0.2,0.5,0.6]

I can expand the equation to write as $y[n] = x[n] + 0.2*x[n-1] + 0.5*x[n-2] + 0.6*x[n-3]$

My problems are:

From a time series or statistics point of view, FIR filter = Moving Average and IIR = Autoregressive model. I cannot understand how to relate the delays, 1,2,3 denoted by the order to the lags terminology often used in statistics. The impulse are known as the coefficients of the model and how to relate the impulse to the coefficients?
Why or rather what is the intuitive meaning of the coefficients becoming zero?
When does that happen/some real life examples of applications?


texts: What is the general feature of a time domain signal that gives a step like spectrum?

texts:
I'm trying to study a nonlinear system by sending a pulsed signal into the system and look at the response spectrum. The signal I send to the system is something like

and the response spectrum from I get from the system is something like (note the spectrum is plotted log scale, ie, (frequency, log of absolute value of Fourier transform at that frequency))

where we can see that there is a stair like structure in the spectrum. So what is the general feature the time domain response signal should have in order to give a stair like spectrum like this? And have you ever encountered a spectrum like this, and is there something in general we can say about the system just from the stair like spectrum?

Could you give a model of this behavior that generate stair like spectrum?


texts: Apply Singular Value Decomposition on Gabor filters

texts:
According to my question in this link, I applied (with a big help from lennon310) the separability of Gabor filters of any orientation.
The only disadvantage of this method is the convolution in the complex domain. So we can treat  the Gabor ﬁlter G(x,y) as a two dimensional matrix that has a low rank, since the ﬁlter is fairly smooth at lower frequencies. Singular value decomposition can then be used in order to decompose the ﬁlter into a linear sum of real separable ﬁlters.

Here, ui, vi are the columns of the orthogonal matrices U,V. s_i is the singular value. The convolution of the ﬁlter G(x, y) with the image I(x, y) can now be approximated using:

Please I need a help about how can I edit my code in this link in order to apply that in MATLAB. Any help will be very appreciated.


texts:
I need to implement the RLS algorithm but it's for matrices instead for vectors, I have made the below code, but still something wrong is not working well,

EDIT:

The code should be done as below,

But how can we process the line of g(:,n)=P*colo*((lamda + colo'*P*colo).^-1); This notation is valid if H is a vector, and y(n) is scalar, but now we have H is matrix and y(n) is a vector.

    clear all; clc;
N = 1000;                    %% Number of samples
H = randn(64,31);           %% The input matrix
c1 = randn(31,1);           %% The outpu desired vector;
p = length(H);
g1 = zeros(p,1);            %% initial RLS weight;
lamda=1;
sigma=1;
P=(sigma^-1)*eye(p);

for n = 1 : N
    y(:,n) = H'*g1;
    alpha(:,n) = c1(:,n) - y(:,n);                     % Instantaneous error of RLS
    g(:,n)=P*colo*((lamda + H'*P*H).^-1);   % Gain vector
    P=(lamda^-1)*P - g(:,n)*H'*(lamda^-1)*P;    % RLS intermediate term
    g1 = g1 + alpha(n)*g(:,n);              % Weight update rule of RLS
end

As you know the normal RLS algorithm is working when the input data, which is H in my above code, is vector. so it will update the loop for as below:

for n = 1 : N
    y(n) = H'*g1;
    alpha(n) = c1(n) - y(n);                     % Instantaneous error of RLS
    g(:,n)=P*colo*((lamda + H'*P*H).^-1);   % Gain vector
    P=(lamda^-1)*P - g(:,n)*H'*(lamda^-1)*P;    % RLS intermediate term
    g1 = g1 + alpha(n)*g(:,n);              % Weight update rule of RLS
end

That will work normal since y(n) and alpha will be scalars in that cases. but for my case, that will be vectors.
I think, maybe we can convert the matrix H into a vector and then when finish the loop we convert it back to matrix ? but I'm not sure and I don't know how to do it, it's just an opinion.

Your help is highly appreciated. I really need to understand that point.

Thank you


texts: Digital low pass filter vs Kalman filter

texts:
I have experience with the design of FIR, IIR digital filters. I also know about the Kalman filter, but I am not skilled at using them. Consider the case of a low frequency signal from discrete samples and the signal is corrupted by high frequency noise. It seems a digital low pass filter and a Kalman filter are two ways of removing the high frequency noise. When is it best to use a digital low pass filter, and when is it best to use a Kalman filter?

* EDIT *
More specifically, it seems a FIR filter with linear phase or an IIR filter with nearly linear phase might be a better estimator than a Kalman filter in some cases. This might be true when the desired signal is low frequency and the noise is limited to the upper frequencies. A Kalman filter is designed for Gaussian noise, and I described a case where a linear phase digital low pass filter would work very well.


texts: how can we calculate Mathieu wavelet filter coefficients from its design equation? a detailed explanation will be highly appreciated

texts: Elliptic-Cylindrical Wavelets: The Mathieu Wavelets


texts:
I am trying to demodulate a FM signal with a RTL-SDR and python using the arctan method. From here I know that
$$
\frac{d}{dt}\arctan(\frac{q(t)}{i(t)}) = \frac{i(t)\frac{d}{dt}q(t)-q(t)\frac{d}{dt}i(t)}{[i(t)]^2}
$$

and so far I have successfully been able to demodulate using the RHS and the following python code

def discrim(x):
    X=np.real(x) # X is the real part of the received signal
    Y=np.imag(x) # Y is the imaginary part of the received signal
    b=np.array([1, -1]) # filter coefficients for discrete derivative
    a=np.array([1, 0]) # filter coefficients for discrete derivative
    derY=signal.lfilter(b,a,Y) # derivative of Y,
    derX=signal.lfilter(b,a,X) # " X,
    disdata=(X*derY-Y*derX)/(X**2+Y**2)
    return disdata

However I would also like to get the same results using the LHS and naively assumed that python code as follows would do that

def discrim_bad(x):
    X=np.real(x) # X is the real part of the received signal
    Y=np.imag(x) # Y is the imaginary part of the received signal
    b=np.array([1, -1]) # filter coefficients for discrete derivative
    a=np.array([1, 0]) # filter coefficients for discrete derivative
    x = np.arctan2(Y,X)
    der=signal.lfilter(b,a,x) # derivative
    return der

Unfortunately this is not the case. I assume it has something to do with the discretization and most likely the discrete derivative approximations but I am not quite able to figure it out.


texts: Why does the Fourier series not include phase information?

texts:
From my understanding, the Fourier series is a way to describe an arbitrary continuous signal in terms of sinusoids of varying frequencies and amplitudes, as shown here in Equation 1.

Why are the phases of the sinusoids not considered?


texts: Formula for designing Lynn's low pass filter

texts:
For the second order Lynn's low pass filter, the general form of the transfer function is:
$$
H(z)=\frac{(1−z^{-m})^2}{ (1−z^{-1})^2 }
$$
where $m$ is a positive integer. The gain for this is $m^2$, time delay is $m-1$ samples and the nominal frequency corresponding to the location of the lowest frequency zero is $f_l = f_s /m$. The frequency response is given by:
$$
|H(e^{j\omega T})|=\frac{\sin^2(m\omega /2f_s)}{\sin^2(\omega /2f_s)}
$$
At the 3 dB cutoff frequency, for unit gain we have:
$$
|H(e^{j\omega_c T})|=\frac{\sin^2(m\omega_c /2f_s)}{m^2\sin^2(\omega_c /2f_s)} \\ \therefore  \frac{1}{\sqrt2}=\frac{\sin^2(m\pi f_c /f_s)}{m^2\sin^2(\pi f_c /f_s)}
$$
Now, to design the filter for our required cutoff frequency $f_c$ and sampling frequency $f_s$, we need to find the value of $m$ corresponding to these inputs.
So how can we get this equation to the form $m = f(f_c, f_s)$, so that we can avoid getting the $m$ value by trial-and-error or numerical methods? Or is there some other quick approximation formula that we can use?


texts:
Here's my basic understanding of how to implement a convolution reverb:

I measure the IR of a space using a sine sweep. The raw mic output is stored as a Voltage-Time data set. Let's assume that this recording is 5s long.

I run a fft on the recorded mic output to derive the Mag and Phase at every frequency bin. Call this the "IR Buffer".

Next, I have my music source. Let's say it's a vocal track that's playing in real-time. I record this vocal track into a 5s buffer.

Similar to step 2, I run a fft to extract the vocal track buffer's Mag and Phase values at every frequency.

Then, I multiply the Mag and Phase values of the vocal track with that of the IR Buffer at every frequency bin.

Finally, I run a ifft on the multiplied buffer and play the output as the reverb'd vocal track.

Is my understanding above correct?
What I'm very uncertain about is that it looks more like I've created a simple EQ rather reverb.
More specifically, I can't figure out how to extract and apply the time-variant Mag and Phase characteristic of the reverb space. For example, there's a certain decay time for 100Hz in the space, which means that Mag and Phase varies with time for this frequency. How is this time-variant component added to the vocal track?


texts: Spatial Aliasing - Wrap Around F-K Spectra

texts:
I've been using F-K Filter, for a while, but I guess I never had good basic understanding about it Math. Someone asked me what is the cause of frequency wrap around in F-K Spectra plot ? I know it's because of aliasing. But if somebody please elaborate more on the cause of this wrap around? Simple Math explanation perhaps.

Thanks

:)


texts:
I came across a concept that says we can design LPH and HPF from APF when the transfer function of the two APFs given as

$$\begin{align}A_0(z)&=1\\
A_1(z)&=\frac{-a+z^{-1}}{1-az^{-1}}\\\\
H_{LP}(z)&=A_0(z)+A_1(z)\\
H_{HP}(z)&=A_0(z)-A_1(z)\end{align}$$

i tried to implement this in matlab but i did not get the desired frequency response..kindly help me with the matlab code for the same..


texts: Coherent Sampling vs Windowing for 50/60 applications

texts: Let's say I need to sample 50/60 Hz signals to find RMS. I can use for example 1kHz sampling frequency. That means 100ms sampling is 5 cycles of 50Hz and 6 cycles for 60 Hz. That's coherent sampling. But does this account for frequency variations? I don't have zero cross detection in hardware, maybe if I do it in software (microcontroller) it may help. I've heard that windowing the signal help reduce discontinuities in case I don't exactly get a integer multiple of the frequencies during sampling. Should I use windowing? What other approaches can be implemented?


texts: Using FFT for implementing a filter bank

texts:
I am trying to use Fast Fourier Transform (FFT) for decomposing an audio signal into 8 sub-bands according to this link but the problem is the frequency response of the result contains only the first sub-band of the signal. Here is my code in Jupyter Notebook:
%matplotlib notebook

from scipy.io import wavfile      # For reading the wavefile into a numpy array
from scipy import signal as st    # Signal Processing toolbox, used for filter design
import matplotlib.pyplot as plt   # Used for genearating the plots
import numpy as np                # Standard vectorized programming library
import scipy.fftpack              # Calculating the FFT
from IPython.display import Audio # Playing audios

fs, data_ = wavfile.read('Track32.wav')
data_ = data_[:,0] # Select the dara of the first channel
print('Sampling frequency of th read .wav file is: {}.'.format(fs))

# Find the next power of two for zero padding
next_pow2 = int(np.ceil(np.log2(len(data_))))

data = np.zeros(2**next_pow2)
data[(2**next_pow2 - len(data_))//2 : (2**next_pow2 + len(data_))//2] = data_

data_16 = np.reshape(data,(len(data)//16,16))
data_subbands = scipy.fftpack.fft(data_16)
data_f = scipy.fftpack.fft(data_subbands,axis=0)

xf = np.linspace(0,fs/2,data_f.shape[0]//2)
fig,ax = plt.subplots()

bands = ['0 - 2 KHz','2 - 4 KHz','4 - 6 KHz','6 - 8 KHz',
        '8 - 10 KHz','10 - 12 KHz','12 - 14 KHz','14 - 16 KHz']

for pli in range(np.shape(data_f)[1]):
    ax.plot(xf,(2/data_f.shape[0])*abs(data_f[:data_f.shape[0]//2,pli]))

ax.set_title('Frequency Response of the Filtered Signal')
ax.set_ylabel('|Y(f)| [dB]')
ax.set_xlabel('Frequency[Hz]')
ax.legend()
plt.show()

Using the general method for filter banks, I came to this conclusion:

But after using this code above, I came to this figure:

Is this the right way to use FFT for implementing a filter bank?


texts: What is the effect of padding zeros to FFT on time domain ? (opposite of what is usually done)

texts: As above. I know that zero padding in time increasaes frequency resolution but I'm having a hard time to understand the opposite. How would you interpret that?


texts: How to undo the quantization step in DCT?

texts:
How can I do point by point multiplication of the quantized DCT coefficients with the quantization matrix before computing the inverse DCT for each block? I want to modify the reconstruction algorithm to approximately undo the quantization step shown here:

for i = 1:Blocks
    for j=1:Blocks
        irange = (i-1)*8+1:(i-1)*8+8;
        jrange = (j-1)*8+1:(j-1)*8+8;
        quantblock = round(y(irange,jrange)./Q);
        quant_DCTblocks(irange,jrange) = quantblock;
    end
end
figure
imagesc(quant_DCTblocks)
colormap(gray)
set(gca,'Xtick',[])
set(gca,'Ytick',[])
axis('square')
title('quantized DCT Coeffs')


texts: noise cancellation of non-repetitive noise

texts:
There are a new noise cancellation headphones that almost fully cancel human noise.

What kinds of algorithms could achieve that ?

EDIT:by human noise , i mean talking. In general, previous audio cancellation headphones didn't work that well for attenuating talking.


texts: Connecting edges detected by an edge detector

texts:
I have a binary image obtained from a canny edge detector. The edges are not nicely detected in the center and I need to join them. The connection of edges is orientation and neighborhood dependent. I need to connect the vertices if they are below a threshold (say less than 5 pixel distance). If my almost linear features are nicely oriented, this threshold should be relaxed. (This is the best case.)

I tried the Hough transform, but it didn't work for me because I don't have straight lines. Dilation and erosion are also not good; they make the images messier.

The approach I am trying is to first detect vertices and nodes (in MATLAB with bwmorph) and then make leaves as a separate feature. This is done as follows:

Scanning in a 3x3 window, look for neighbors.
Traverse the complete connected object.
Try to fit a line (or maybe a quadratic polynomial).
Check feature-by-feature if it's worth a connection or not.

The implementation is not straightforward since the decision-making part where vertices need to be connected is tricky.


texts: Sampling Theorem: How to know the value between two samples of a Signal

texts:
According to Sampling theorem, in order to reconstruct a signal we need to sample it at the rate => twice the highest frequency component of that signal. (provided signal is band limited).

Let's say, we have a signal with f= 2MHz (highest freq component), so we will be sampling it at 4MHz or more as Sampling theorem say's.. and we will have N Samples in signal.

Now, what if I want to know the exact value of a sample (information) between two samples..??

How can we find the value at every possible instantaneous.. in a signal..when we know we only have N samples after sampling..??

I know we can increasing sampling rate to have more sample. But is there any other way to do it..?

This question was asked to me in an interview at NCRA-TIFR, Pune, India.


texts: Which books do you recommend that provide comprehensive overview of delay compensation techniques in control systems?


texts: Best Approach on How To Sample and Process

texts:
I'm looking for some advice, I am hoping one of the many forum experts could help me with some advice on the sampling and processing of signals. I am using an STM32F4 processor with 12 bit ADC and DSP core. Presently, I can use either CMSIS DSP or Siglib libraries. Here's a point description of the problem space:

I can sample at whatever rate is required: in front of the ADC is a 10th order analog amplifier and LPF, Fc = 55Hz. Amplifier can be set to provide ideal ADC range of 3Vpp

Sinusoidal incoming signals are in the range of 10Hz to 40Hz that can be set for any duration in time. Ideally, in the range of 200ms to 300ms (can be changed)

In between these tones, there is dead space of (ideally) 50ms to 100ms (can be changed)

There would be a series of these tones with a total time of about 4, or possibly 5 seconds (can be processed after last signal arrives, ideally real-time)

I need the best amplitude and frequency accuracy possible, ideally about 2dB and 0.5Hz. The environment is noisy, but incoming signals can be kept 20dB over noise.

I am not looking for a code handout. Just a best approach algorithm and ideally some specifics such as sampling rate (is decimation required), buffer size, to use overlap-add/sliding window, signal processing chain - FIR required (how many taps), FFT method/size?

It may not be complex to many of you, but I haven't finished all the theory I need to try an implementation without spending a long time,(I have all the books but there's a lot of reading yet!) and the problem requires me to try for a solution now.

Thank you for your help!
Roger


texts: Converting a parallel filterbank to cascade structure

texts:
I have a bank of 12 bandpass filters (biquads) that I use for spectral modelling.

Each filter is centered around a specific "mode" frequency, and has its own resonance.

Until now I have implemeted this filterbank in parallel structure, i.e. by summing the output of each filter aplied to the same input. I also normalize the gain of the filters to restrain the output range.

However, to make proper use of the Cortex M4 CMSIS-DSP functions and thus improve (I hope) the performance of my code, I would like to try and make the filterbank in a cascaded structure. I have tried on Matlab to directly change the structure without changing the coefficients, but the sound is not the same at all and the output range isn't normalized at all.

Is there an easy way of converting the coefficients of each biquad from the parallel form to the cascade form?

Additional info :

Here is how I calculate my coefficients :

$$
\begin{split}
ff0 &= \sqrt{\frac{1-R^2}{2}}\\
ff1 &= 0\\
ff2 &= -\sqrt{\frac{1-R^2}{2}}\\
\\
fb1 &= -2\times R \times cos(\omega) \\
fb2 &= R^2\\
\end{split}
$$

Where $R$ is the resonance of my filters and $\omega$ is the frequency-normalized pseudo-pulsation.

In Matlab, the main loop for filtering the input is :

for i = 1:length(fc)
    out = out + filter(ff(:,i), fb(:,i), in); % Parallel
    outCasc = filter(ff(:,i), fb(:,i), outCasc); % Cascade
end


texts: Why doesn't law of large numbers apply to this stationary time-series?

texts:
There's a paragraph in Wikipedia that states the following:

  Let Y be any scalar random variable, and define a time-series $\{X_t\}$, by
  $$X_{t}=Y\qquad {\text{ for all }}t$$
  Then $\{X_t\}$ is a stationary time series, for which realisations consist of a series of constant values, with a different constant value for each realisation. A law of large numbers does not apply on this case, as the limiting value of an average from a single realisation takes the random value determined by $Y$, rather than taking the expected value of $Y$.

I cannot understand the reason for the law of large numbers not being applicable. Any explanation would help.


texts: Sinusoidal liftering in implementations of MFCC

texts:
Some implementations of MFCC apply sinusoidal liftering as the final step in calculations of MFCC. It is claimed that speech recognition can be significantly improved. For instance, if $\text{MFCC}_i$ is a cepstral coefficient, and $w$ is a lifter, then
$$\widehat{\text{MFCC}_i}=w_i\text{MFCC}_i$$

is a liftered cepstral coefficient, where $w_i$ for sinusoidal liftering is defined as:

$$w_i=1+\frac{D}{2}\sin\Big(\frac{\pi i}{D}\Big)$$

When I look at the equation, I understand the sinusoidal function has a shape such that its maximum is in the middle and approaches to zero at edges. Therefore, cepstral vector's first and the last coefficients are is reduced to zero while the middle one is intact.

Why is liftering applied and how does it improve the speech recognition?


texts: Subtract background from 2D signal using spectral methods

texts:
I have a 2D signal, with intensity as a function of two angles (alpha and beta), as shown below.

This contains a background, which makes the signal's base-surface locally convex-concave. I would like to identify and subtract the background such that the bottom of the signal becomes approximately planar and further processing can be done.
My initial thoughts were to fit row-wise/column-wise baselines and individually subtract the background from each line. However, I think signal processing methods may be better suited for this task. I have calculated the FFT for the signal using SciPy/Python, but don't find any obvious frequencies that I should remove.
I'm not even very clear about whether I should be looking to remove the highest frequencies, the lowest frequencies, or both. What would be a good way to proceed from here, using Python? Is the power spectral density better to look at rather than the FFT itself?
EDIT :
The signal comes from a diffraction experiment, where there is a source, a specimen and a detector. The specimen is tilted from 15 degree to 90 degree in steps of 5 degree. At each tilt, it is also rotated through 360 degrees in steps of 5 degree. The tilt and rotation are the angles alpha and beta shown in the figure. Corresponding to each alpha/beta combination, the detector receives some signal, which we are plotting as intensity.