texts: Auto/cross correlation of a cosine in the context of an adaptive beamformer

texts:
I am currently reading chapter 13 of 'Adaptive Signal Processing' by Widrow and Stearns, so if anyone happens to have a copy of the book to hand it could be helpful.

I am reading the chapter on "Introduction to Adaptive Arrays and Adaptive Beamforming" and am struggling understanding the auto-correlation/cross correlation between the primary and reference signals at a delay of one (as shown in the diagram below) and in the more general case (where the signal is arriving at an angle), where the general cross correlation vector is described by:

$$
\mathbf{P} = \begin{bmatrix}
\phi_{dx}(0) \\ \phi_{dx}(1)
\end{bmatrix} = \begin{bmatrix}
E [C \cos k\omega_0 \cdot C \cos (k + \delta_0)\omega_0] \\ E [C \cos k\omega_0 \cdot C \sin (k + \delta_0)\omega_0]
\end{bmatrix}
$$

where $k$ is the sample number, $C$ is a constant amplitude, $\delta_0$ is the time delay of arrival between the primary and reference sensor and $\omega_0$ is the angular frequency.

The book assumes that the refence signal to the array is a cosine of the form:
$$
\text{reference signal} = C \cos [(k + \delta_0) \omega_0].
$$
(The above equation is 13.3 in the book).

My question is, I don't understand why, for $\phi_{dx} (1) $ the 'lag' part of the equation becomes a sine rather than a cosine with a shift 1 (sample)? It's as if there is a shift of $-\pi/2$, hence the sine, but I can't see where this is from? Why for a lag of 1 does this cosine seemingly change to a sine?

My understanding of cross (or auto) correlation is that for $\phi_{dx} (1)$ you compare one signal with another with a shift of 1 sample (hence a lag of 1).

I understand that a sine is the quadrature of a cosine and understand why these equal zero (as the book describes; "$\phi_{dx} (1) $, is zero because it represents the correlation of (13.3) with its quadrature (sine) component"). My question is why is a sine introduced for a lag of 1?

Diagram for completeness:

Many thanks


texts: How to, efficiently, identify telephone number from the key pressed data audio file?

texts:
A very common homework question. You are given an audio file(.wav) which has the recording of the keys pressed during a telephone call. You have to identify the telephone number from this audio file when you are given the sampling frequency & the duration for which one key is pressed.

From the given sampling frequency & duration of key press, I could manage to write a MATLAB code to identify the telephone number by taking (samp_freq * key_press_time) number of samples at a time, doing fft on it & finding the bin that gives me a peak. I used WAVREAD command to read the audio file.

One problem I noticed is that, if the duration of key press is different(1st key pressed for 1 sec, 2nd key for 2 sec, 3rd key for 0.5 sec & so on), my technique of separating the data & finding peak will fail.

How to get around this problem of random duration key press? I could think of 2 inefficient ways to do this.

During 1st iteration, take entire data & plot fft. If (no. of peaks > 1), take half data at a time plot fft. Keep on doing this until there is only 1 peak in the plot. But this technique will take a lot of time
WAVREAD the file, each key press data is separated by 0's. Thus take data between the 0's one at a time, plot fft & find peak & frequency. But this looks too simple to be used.

So is there any other effiecient method of finding the telephone number when the key press duration is random?


texts: Deconvolution to Remove Gaussian Blur in 1D Signal (Wiener Filtering?)

texts:
I've got a set of biology data that I'm trying to denoise (effectively, a population statistic can only be measured convolved with a gaussian of known width)
My problem is this: I can measure (f*g), and I can measure g(x), but I need to know f(x). Because this is real data with slight blips, I can't just use a deconvolve() function.
I've been told by an EE friend that I should be using 1-D Wiener Deconvolution. The problem is, I have no earthly idea how to implement this. Supposedly, there's a Matlab package for it, but I don't have access to Matlab; and all the Python examples I've been able to find (e.g. scikit-learn) are 2D only.
Is there something I'm missing? Is there an easy formula I can implement?
Please, don't just point me to the Wikipedia page/a textbook for 1D Wiener filtering. I've tried to read several of these, but they look like incomprehensible gibberish (I'm a biologist).
EDIT: Also what's up with point spread functions? I assume my PSF should be the known Gaussian g(x), and yet every example I can find of Gaussian deconvolution with Wiener filtering (all in 2d) use a uniform distribution as their PSF instead of a Gaussian distribution.
Thanks!


texts: Differences between FFT from acceleration, velocity and displacement data of a dynamic system

texts:
I am new to vibration mechanics and signal processing. Recently, I’ve performed an experimental test with an accelerometer mounted on a plate (shaker). The plate vibrates with exact frequencies 10 and 20 Hz. The purpose of the experiment is to extract the acceleration vs. time and FFT diagrams. I have a basic and fundamental question.

What are the differences between applying FFT to acceleration, velocity or displacement diagram of a dynamic system? Does FFT give the same peaks for the three mentioned diagrams so that I assume that the peaks are natural frequencies of the vibrating system (10 and 20 Hz in this problem)?

I would be appreciated if you could help me with these questions.


texts: Detect quasi-parallel lines with spurs in a binary image

texts:
I'm looking for a way to remove quasi-parallel curves from the net, as depicted in the figure. Not only the original nets (upper row in figure) contain many spurs, but also the curves to be detected are not necessary straight lines. Ideally, I'd like to keep only the border of each such a region (lower row in figure).

I've came across solutions like hit-or-miss with line structuring element, Hough transform and Radon transform - but they were aimed at either detecting long straight lines or arbitrary shapes (like ellipses).

My idea was to represent a net as a graph, linearly approximate each edge (and split it into smaller segments if it deviates from the approximated line too much), for each pixel store the slope of an approximated edge in a matrix, and finally for each pixel I'd search such a matrix in a given neighborhood for other pixels (belonging to a different edge) with a similar slope.

Could you recommend me any method already used to accomplish such a task or tell me whether my approach seems reasonable?


texts: Finding Fourier transform of a discrete signal from its Z-transform

texts: Is it possible to find the Fourier transform of a discrete signal if you know its $\mathcal{Z}$-transform of?


texts: How does a sound card translate amplitude values to sound when playing WAV/PCM/FLAC files?

texts: I'm reading the documentation for WAV files (and importing them in C++)It seems like 44100 amplitude (this can change, and is specified by the WAV header) values get sent to the sound card every second. What I don't understand is that these amplitudes can be negative. I'm confused as to what these amplitude values represent. How can they be negative and positive? What happens when you send an amplitude value to a sound card to produce a sound that you can hear?


texts: I have complex baseband signal and its output when passed through power amplifier. How to plot the AM-AM and AM-PM plots?


texts: I am trying to convert a signal into another unit, but the error keeps increasing during each iteration

texts:
I have a signal, A, and I want to convert it into another unit by using a formula:

$$s_i = \frac{(1+x_i)}{(1-x_i)}\cdot s_{i-1}$$

where $x_i$ is each sample of the signal. I am using a loop where the last value of $s$ is stored and used for the next iteration (as the new value of $P$). It goes on and on until the last sample of the signal.

Signal B is the result of the calculation. However, I noticed that the error keeps increasing during the calculation. As shown in the figure, the red line is how I want the signal to be.


texts: How to Classify a Kernel as Low Pass Filter (LPF) or High Pass Filter (HPF)? How to Transform an LPF Kernel into HPF Kernel?

texts:
The filter impulse response is [1/9 1/9 1/9;1/9 1/9 1/9;1/9 1/9 1/9].
how do i determine whether it is an Low pass filter or a high pass filter. If this is a low pass filter how do i convert it to a high pass filter or the opposite.


texts:
My book give me two signals to demonstrate that the temporal translation does not alter the energy and area. It gave me

$$ x(t)=\operatorname{sinc}(t) $$

and

$$ s(t)=x(t-T)$$

and I found that energy ( with rayleigh theorem ) and area are 1 for $x(t)$. For $s(t)$ I found an area=1 but now , with Rayleigh I should demonstrate that.

$$  \int\limits_{-\infty}^{+\infty} | \operatorname{rect}(f) \,e^{- i 2 \pi f T } |^{2} \ \mathrm{d}f  $$

As I made for $x(t)$,

$$ \operatorname{rect}(t) = |\operatorname{rect}(t)|^{2} $$

but solving the last integral this didn’t give me 1. Thank you so much for the help.


texts: Cross correlation with short or non periodic signals?

texts:
I am a total newbe to signal processing and i have a problem which may be  quite trivial to most of you guys.

I want to align two short sensor signals (speaking of max. 1000 samples recorded with 20Hz). The signals are sampled with the exact same sample rate and are not necessarily periodic. I use octave to make first experiments and i recognized some confusing behaviour of the xcorr() function in the 'signal' package.
If i correlate two short signals the resulting lag value is not really correct. The longer one of the signals is, the better the lag value calculated by the xcorr() function.

Here is what i tested:

pkg load signal

clear all

%generate two short signals
x = -pi:0.05:pi;

sigA = sin(x);
sigB = cos(x);

%plot the original signals
figure 1
subplot(2,1,1);
plot((0:length(sigA)-1), sigA, 'r');
title("sin(x)");
xlabel("Sample index");
grid on;
subplot(2,1,2);
plot((0:length(sigB)-1), sigB, 'b');
title("cos(x)");
xlabel("Sample index");
grid on;

%calculate the cross correlation between signal B and signal A
[c, lag] = xcorr(sigB, sigA);
[~,i] = max(abs(c));
shift = lag(i) %the shifting offset between the two signals

%align the signals
sigA = sigA(abs(shift)+1:end);

maxLength = max(length(sigA), length(sigB)); %for sample index axis range

%plot the aligned signals
figure 2
subplot(2,1,1);
plot((0:length(sigA)-1), sigA, 'r');
title("sin(x) aligned");
xlabel("Sample index");
axis([0, maxLength, -1, 1]);
grid on;
subplot(2,1,2);
plot((0:length(sigB)-1), sigB, 'b');
title("cos(x) aligned");
xlabel("Sample index");
axis([0, maxLength, -1, 1]);
grid on;

The code produces the following figures and calculates a shift of -27 samples:

As you can clearly see, the signal A (sin(x)) is not perfectly aligned!
I would have to change the alignment call to

%align the signals
sigA = sigA(abs(shift)+6:end); %where does this offset (1+5) come from?

to align them perfectly:

Now if i generate a longer signal for sigA, this offset disappears completely:

x = -100:0.1:100;
y = -pi:0.05:pi;

sigA = sin(x);
sigB = cos(y);

Then the situation looks like this:

As you can see in the image above, the two signals are perfectly aligned without shifting signal A manually as i had to with the shorter version of signal A.

Now here is an actual signal from my sensor (where signal A is a stored signal and signal B is a subset of signal A, the x axis shows the time in ms):

Same Problem appears here. I would have to shift signal A by 29 samples to the right to align them perfectly.

Why is that? Is cross correlation not the right thing to do for short or non periodic finite signals? Or is there a way to calculate the additional offset to the signal  shift?


texts: Difference: Ergodic channel capacity and spectral efficiency

texts: Can someone explain what is the difference between the capacity and ergodic capacity? Ergodic capacity and spectral efficiency?


texts:
In terms of signal analysis (especially frequency decomposition) would it be bad to use cyclic autocorrelation (AC) instead of acyclic/linear AC?

My use case:
I'm trying to characterize events in a sensor stream by periodicity/frequency. At the moment, I'm using DFT spectrum (spectral density) and AC of sliding windows.
On the one handside it's algorithmically way cheaper to calculate the cyclic AC via FFT than acyclic (pad with zeroes and do double length FFT).
However, the cyclic AC looks similar (though not equal) to the acyclic for the first N/2 lags and is self-symmetric - thus at most half the information of acyclic AC can be contained.

My motivation:
If the cyclic AC would be sufficient for detection, that would mean better performance compared to linear AC (in terms of memory usage and computation time). I then could also calculate AC inplace (better parallelization potential).

Rephrased question:
Is there any rule when cyclic AC should NOT be used?
Or perhaps, are there known cases, when cyclic AC might actually be better than linear AC?

Thanks in advance :-)


texts: Different results for different orders of estimating AR model using Yule-Walker equations

texts:
I'm trying to use MATLAB to estimate the AR parameters to the following filter:

$$H(z) = \frac{1}{1-0.5z^{-1}+0.25z^{-2} -0.25z^{-4}}$$
As I can see, the process at the output of this filter depends on the previous 4 samples. So, i figured that i would need 4th order model of Yule-Walker equations.

I run a MATLAB code solving Yule-Walker equations for 4,8,16,64 order and when I compared the 4 options, I got that the higher the order, the higher the MSE error. How can i explain this fact?

Here is a plot I got when compared $\lvert S(w)-S_{\rm estimated}(w)\rvert^{2}$.


texts:
In the scale-space theory the scale-space representation of the signal $f(x), x = (x_1, ..., x_d)$, (in case of image $d = 2$) is given as:
$L(x, y; t) = g(x, y; t) * f(x, y)$ where $g(x, y; t)$ is a gaussian kernel with parameter $t$ and $*$ is a convolution. By changing the $t$ parameter we receive a more or less smoothed image. As the result coarser representation (parameter $t$) will not contain small objects or noise.

The main point is to find a way of scale-invariant feature detection, right? So that for some image at it's reduced in size copy the features like keypoints will be detected correctly, even if size is different, without finding other noise keypoints.

In the paper they are using the $\gamma$-normalized derivatives.
$\delta_{\xi, \gamma-norm} = t^{\gamma / 2} \delta_x$.
What is the meaning of using the $\gamma$-normalized derivative, how does it helps in scale-invariancy?
From this image we can see that at near the same positions the different keypoints found (different in size). How is that possible?

If you can explain the step-by-step algorithm of scale-invariant feature detection, this would be great. What is actually done? The derivatives can be taken by $x, y$ or $t$. Blob can be detected by taking the the derivative of $L$ by $(x, y)$ variables. How is the derivative by $t$ is helping here?

The paper I was reading is: Feature detection with automatic scale selection


texts: Equalizer coefficients and channel coefficients

texts: Blind channel equalization methods equalize the channel without using the source data and without knowing the impulse response of the channel. Consider the channel to be a single input single output FIR system. The popular methods for blind channelequlization are Constant Modulus Algorithm and Least Mean Squares. Certainly, the equalizer coefficients are the not the estimated channel impulse response. It is not clear to me that once the equalizer weights are estimated, how does one get the estimates of the channel coefficients/ impulse response? I am missing a fundamental step somewhere, can somebody please help how we use the equalizer weights to estimate the channel?


texts: Filter design for an unusual EEG experiment

texts:
I'm designing a quite unusual (in terms of signal analysis) EEG experiment, which forces me to design my own filter. I have never done this before, so kindly please check my work and suggest corrections.
The analysis consists of cutting a narrow frequency band (2Hz) from 6-second EEG recordings(4*) (sampling frequency 256 Hz). I expect that after such filtering, I will receive a "somehow oscillated" signal (not in a strict mathematical sense, just with notable local maxima (peaks) and local minima (troughs))(5*). Now I'm interested just at the very end of this 6-sec window, the last dozen or so data points. I want to determine (by the simple algorithm provided by myself) if the signal ends "near the peak" or "near the trough" or in either of the above. So, all 6-sec windows are compressed to one of the three categories. You can think of this as a very rough estimate of the phase at the end of the window.
The tricky part is that I need to analyze overlapping frequency bands ( 4-6 Hz, 4.25 - 6.25 Hz, 4.5 - 6.5 Hz, etc.), so it will be great to have an extremely steep roll-off even at the price of amplitude disturbance (most of this information will be lost anyway) - that's why I came up with an elliptic IIR filter. Zero-phase shifting is also crucial so, forward-backward (filtfilt) filtering is applied. However, the devil is in the details (filter order, maximum ripple, minimum attenuation), and I must honestly admit that I determined these parameters by trial and error.

from scipy import signal
from mne.viz import plot_filter

sfreq = 256.
f_s = 6.
f_p = 8.
flim = (f_s - 1.0, f_p + 1.0)  # limits for plotting

nyq = sfreq / 2.  # the Nyquist frequency is half our sample rate
freq = [0, f_s, f_s, f_p, f_p, nyq]
gain = [0, 0, 1, 1, 0, 0]

ftype = "ellip"
order = 8
sos = signal.iirfilter(order, [f_s / nyq, f_p/nyq],
                       btype='bandpass', rp=5, rs=35, ftype=ftype, output='sos')
plot_filter(dict(sos=sos), sfreq, freq, gain,
            f'{ftype} order={order}', flim=flim, compensate=True)

Filters Visualisations:

So my questions are as follow:

Is my filter "reasonable"?
Can I do better?
Can you recommend to me some learning materials for a filter design.
(*) This is obviously cut from continuous EEG signal - can it be just cut (boxcar window), or should I apply something different (eg. Hamming window)?
(*) Is my expectation justified?


texts: Tips for programming scientific articles

texts:
I have to reproduce simulations of scientific articles in Matlab.

Do you have any tips for programming scientific articles? What are your methods when you have to reproduce simulations of scientific articles?


texts:
I have the following vector:

%create vector of sampling
t=0:pi/100:2*pi;
x=sin(t);

What is the sampling frequency (of this vector) ?

Thanks


texts: Rejecting a particular frequency band using DWT

texts:
My problem is as follows: I decompose the 4th level detail coefficient of an ECG signal further up to the fourth level and I want to reject the approx. coeff. at this level, which corresponds to the powerline interference. I reject this by using upcoef() with 'd' and then specifying only the detail coefficient. Then I want to keep on using idwt() to get back at the previous level , but this is not being possible because of the length mismatch of the new reconstucted detail coeffs. and the original approx. coeffs.
I would like to know if at all it is possible (in MATLAB) to reject a frequency band corresponding to the further decomposition of the detail coefficient, because when you go back up, the reconstructed coefficient lengths are different because of the inherent convolution operation.

I am posting my code below:

[C1, L1]=wavedec(detrendedecg, 4, 'db8');
[cd4]=detcoef(C1,L1,[4]);

[aa1,dd1]=dwt(cd4, 'db8');
[aa2,dd2]=dwt(dd1, 'db8');
[aa3,dd3]=dwt(dd2, 'db8');
[aa4,dd4]=dwt(dd3, 'db8');
newdd3=idwt([],dd4,'db8');

newdd2=idwt([aa3],[newdd3],'db8');
newdd1=idwt([aa2],[newdd2],'db8');
newcd4=idwt([aa1],[newdd1],'db8');

the original signal is detrendedecg. It is decomposed upto level 4 to get cd4 as detail coefficient. cd4 is decomposed again up to 4 levels to obtain dd4  as the detail coefficient at level 4. I want to reject the approx coefficient at this level, i.e., aa4. But this code doesn't help. It gives a matrix dimension mismatch error on the third line from thee end. What do I do now?


texts: How to Classify Image Filter as Low Pass Filter (LPF) or High Pass Filter (HPF)?

texts:
I am currently reading on the topics of image filtering and convolution kernels. While the process of convolution is clear to me, I do not understand, how one arrives  at the statement, that a certain filter is a low or a high pass?

How can one derive that? For example, if one averages the value of a pixel respecting its surrounding pixels, does that for all pixels, then each pixel would have a value between the minimum and the maximum value of the previous image. It treads low values and high values equally, however literate often states, that it is a low pass filter.

For that specific case I found this question here. Reading it, however didn't help me, as I seem to lack something to see how these arguments are correct and intuitive. Thus I hope, that a more mathematical approach (can be a general one), would benefit me.


texts:
Question is to find DFT, $G_7(k)$ of $g_7(n)$ = $x(2n)$.
For purpose of illustration only, we can assume that $x(n)$ can be represented by envelope shown in corresponding figure.

Shown in Figure Below, is solution to this problem.

But, I am unable to understand this Solution.
If anyone can help me understand this solution, I will be very thankful.


texts: Is there an equation for efficiency of a 2-radix FFT as you raise $k$?

texts: I know the efficiency for an $N$ point 2-radix FFT is $N\log_2(N)$ but assuming $k\leq N$, what if you were looking for the efficiency of calculating $k$ positions of the FFT? Would the efficiency be $k\log_2(N)$?


texts:
I have to quantify the broadband noise of an audio track.

I've seen lots of things around and the only viable option to estimate the broadband noise (without having any information about it) is something like that:

noise=0;
sig=0;
for i=in:fin
   nois(j)=(abs((min(B(:,i)).^2)));
   sign(j)=(abs((max(B(:,i)).^2)));
   j=j+1;
end
snr=var(sign)/var(noise)

Where B is the matrix of the spectrogram.

So, I'm doing FFT analysis of the audio, finding the minimum amplitude for each bin for an approximation of the noise floor and the maximum amplitude for an approximation of the signal. Then I'm doing the variances of both vectors and finally making a simple signal to noise calculation.

I know it's not a very sophisticated method but things previously tried have not worked. This gives me more solid results but it doesn't work for all kinds of genres.

Has someone something more to suggest me to improve this algorithm? Other approaches that work better?


texts: What Is the Best Color Space to Work with for Hand Detection: RGB, HSV, or YCrCb?

texts: I'm working on hand gesture recognition using C#. I'm using OPENCV for skin detection but its performance is not optimal. So I decided to do some filtering based on color space, but I do not know what color space is the most suitable for this purpose?


texts: How can I create a jamming signal for desired frequency?

texts:
I am working on GPS interferences and countermeasures. I came across jamming in GPS receivers.

So I have a  input signal with frequency of 1200 MHz

How can I create a jamming signal that can mask the input frequency. I am working in Matlab.

Matlab has in built function called barrage jammer, which i didn't understand.

Thanks


texts: Solving for a convolved function inside an integral

texts:
I have an imaging system that produces a signal in a pixel, in part, by responding to incident radiance $L(\lambda)$ that passes through a filter with transmission $T(\lambda)$. Lumping all the other parameters under a generic constant I can write the signal:
$S = k\cdot\int L(\lambda)\ast T(\lambda)d\lambda\,$.
I'm trying to find a way to determine, or even estimate $L(\lambda)$. The system is such that I know $S, k$, and $T(\lambda)$. Is there a way to recover the incident radiance? Note, I am not necessarily looking for a closed form solution, approximation is fine.


texts:
I'm an absolute beginner.
I'm doing a project that involves estimating the time delays between the sound signal arriving at different microphones in an array (the project is a bit complicated but I'll find my way through). I'm planning to use a panda board and connect the microphone array to it. But before going about implementing the array, I thought I should first learn how to do some DSP on an analog input. Basically, how would I connect ONE microphone to the panda board and do, for eg, an FFT on the signal?

The panda board has an audio in port but that wouldn't help much when I have to use multiple microphones. Can I connect the condenser mic (the really cheap one ) using 2 wires to a general input pin? How would I read input from it? Will I be able to record sound using such a setup? And I understand that I have to do the Analog to digital conversion (In fact I WANT to do it by myself, so that I would better understand how stuff work).

And now the million dollar question : If I connect a few more microphones to different input pins of the panda board and give them the same sampling rate (I have no idea how but I'll learn it) will that make a microphone array?


texts:
I try to obtain fundamental frequency from cepstrum. I decide to use peak picking algorithm, which pick maximum value of the Cepstrum (c) with sampling rate fs.

function [f0] = spPitchCepstrum(c, fs)
 % search for maximum  between 2ms (=500Hz) and 20ms (=50Hz)
 ms2=floor(fs*0.002); % 2ms
 ms20=floor(fs*0.02); % 20ms
 [maxi,idx]=max(abs(c(ms2:ms20)));
 f0 = fs/(ms2+idx-1);
end

It works but I can not understand why f0 is calculated like:

f0 = fs/(ms2+idx-1);

where in denominator is value consist indexes.

Can you somebody explain me why it is calculate in this way?


texts: relation between luminance and eigen values of an image

texts:
There are some research papers in the literature which apply Singular Value Decomposition (SVD) on an image to decompose the image matrix into U, S and V components. The eigenvalues in the singular matrix S are multiplied by a suitable factor and the image is then reconstructed. The resulting image's luminance is enhanced. To refer to one such paper, click here.

In curiosity, I was trying to study the variation of eigenvalues with respect to the luminance of an image. For calculating luminance, I am using log-average luminance described here. I will call log-average luminance  Key Value from here on.

Experiment:
I began my experiment by building a data set first. I took back-to-back multiple images of the same scene with different exposures.

For each image in such a set, I calculated luminance and obtained eigenvalues for each image. Following is the kind of graph I got for R, G and B channels respectively:

To clear any ambiguity about what these graphs represent, I'll explain briefly here. These graphs are for one set of images (11 images). For each image in the set, I calculated key-value. Next, I calculated average eigenvalues for each channel.
For calculating average I used only the first 20 eigenvalues.

First, I observed one thing about eigenvalues which I don't understand. Not all the eigenvalues increase with the key-value of image. When key-value was more than ~0.85, out of 20 eigenvalues, only the first two or three tend to increase; the rest started decreasing. So this led to a turning point in the graph for images with key-value more than 0.85, so I neglected those images. Does this mean only the first few eigenvalues carry the maximum information about image brightness? Is there any reason for such observation?

Lastly, if I fit a curve to the obtained data points, can I generalize a formula for all images? So that for a given key-value I can get a fair estimate of eigenvalue (maximum or average). Or is there any other recommendation to fit a curve for getting a generalized formula?


texts: Which is correct representation for function of variable?

texts:
I have a function as
$$E=\int_\Omega -\log\big( p_i(x)\big) dx$$
where $p_i(x)$ is density distribution which estimated by Parzen window method.
$$p_i(x)=\frac{1}{\Omega_i} \int_{\Omega_i}K_\sigma\big(I(x)-I(y)\big) dy$$
where $\Omega_i$ is the image domain, with
$$K_\sigma=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(\frac{-z^2}{2\sigma^2}\right)$$
and $\sigma$ is scalar parameter. $I$ is an image $I(x):\Omega \to \mathbb{R}$, $i$ is given.

So, which ones is correct for left side: $E(\sigma),E(x),E(x,\sigma),E(p_i)$ or other form?

Update: I upload the original equation for better visualization


texts: Initial rest condition for the linear constant-coefficient differential equations

texts:
Suppose that system has the input/output relation as follows

$$\sum_{k=0}^{N}a_k \frac{d^ky(t)}{dt^k} = \sum_{k=0}^{M}b_k \frac{d^kx(t)}{dt^k}$$
Where $a_k, b_k \in \mathbb{R}$.

Obviously we need additional conditions to determine completely the input-output relationship for the system and those auxiliary conditions specify whether the system is LTI and causal or not.
The definition of initial rest condition for a system in general is that if $x(t) = 0$ for $t\lt t_0$ then $y(t) = 0$ for $t\lt t_0$.
First of all is it true that the initial rest condition for the all systems implies $y(t_0) = 0$? Also I have seen expressions like $y({t_0}^+)$ and $y({t_0}^-)$ but I don't know how they relate to $y(t_0)$. Another question is that whether the initial rest condition is necessary and sufficient for the mentioned system to be causal and LTI? What's the rigorous proof, then? I searched on the Internet but didn't found a proof or a clear explanation.


texts: Problem with natural logarithmic function in fixed point notation

texts:
I am facing a problem with natural logarithmic function in fixed-point notation.
Let's say
$$x = 0.54, \qquad\text{then}\qquad \ln(x) = -0.616186, \qquad\text{a negative number}  \tag{1}$$
Then in fixed-point notation, the value of $x$ is calculated by (assuming $x$ is $\rm Q8.24$ format)
temp_var = 0.54 * 16777216 = 9059696 (fixed point value in Q8.24 format)

So,
$$x = \ln(\textrm{temp_var}) = 16.01\qquad\text{since}\qquad 2^{24} = 16777216$$
And the float (true / unscaled form) value corrosponds to
$$16.01 / 16777216 = 9.54e-7\qquad\text{a positive value}\tag{2}$$
As we can see, both $(1)$ and $(2)$ are different in both magnitude and sign.
So, I just want to know how to solve these kinds of problems at a fixed point?


texts: Will noise cancellation with playstation eye improve my spectral analysis?

texts:
Im doing some beginners spectral analysis for a homemade reactive RGB LED project and was thinking of ways to improve the accuracy of the audio analysis. Im using a Playstation eye with a Raspberry Pi Zero W (limited performance) and ive read several articles to try to better my understanding. According to this article, the playstation eye has 4 channels, where the 2nd and 3rd are reversed for noise cancellation.
I assume that the 2nd channel is a reverse of the 1st and the 3rd of the 4th. Am i correct to assume that the way to go about implementing these extra channels would be to average the two sums of 1+2 and 3+4 ie. voltage = (c1+c2+c3+c4)/2?

Below is a script i use for testing, only the basics. Note that i use a sampling rate of 20khz instead of 44.1khz due to the limited cpu resources of the pi zero. I also use int16 as wiki states that the ps eye uses 16bit ints. Am i wrong not to use the default float?

The reason i want to improve it is that i havent recorded any frequencies over 7000, which seems a bit odd to me given the vast variety of music ive tested with. Im wondering if the playstation eye is not suitable for such precise analysis, or if my methods are wrong. Is it possible to improve this with noise cancellation, or will it worsen my results? Should i rather just use a single channel or sum channel 1+4?

#!/usr/bin/env python3

import sounddevice as sd
import numpy as np

RATE = 20000
lowest = 200
highest = 0
h_m = 0

def print_sound(indata, frames, time, status):
    global lowest, highest, h_m
    if status:
        print(status)
    if frames == 0:
        return
    data = np.sum(indata, axis=1)
   # data = data * np.hanning(len(data)) # smooth the FFT by windowing data
    fft = abs(np.fft.rfft(data)) # calculate magnitudes from positive fft coefficients
    freq = np.fft.rfftfreq(frames, 1.0/RATE) # calculate positive frequency bins
    max = np.max(fft) # highest magnitude
    freqPeak = freq[np.where(fft==max)[0][0]] # find frequency with highest magnitude
    if freqPeak < lowest:
        lowest = freqPeak
        print("lowest: %d"%freqPeak)
    if freqPeak > highest:
        highest = freqPeak
        print("highest: %d"%freqPeak)
    if max > h_m:
        h_m = max
        print(max)

with sd.InputStream(dtype='int16', channels=4, samplerate=RATE, callback=print_sound):
    while True:
        response = input()
        if response in ('', 'q', 'Q'):
            break


texts:
How to obtain the correct amplitude after the numerical Fourier transform of a signal?
Example: consider an exponential decaying wave $y(x)=e^{-x}\sin(100\pi x)$ with Fourier transform $y_f(x_f)$ which should be a peak with some broadening. I am searching for the appropriate approach to numerically approximate the continuous function $y_f(x_f)$ as close as possible. The calculation of the approximation of the function $y_f(x_f)$ should thus be independent of the number of samples and the range that y(x_i) is defined.
However, when I calculated the Fourier transform for several x-ranges (N in the code), I obtain different results:
fig, ax = plt.subplots()
N=[600,1200,2400,4800]
for i in N:
    N = i
    # sample spacing
    T = 1.0 / 800.0
    x = np.linspace(0.0, N*T, N)
    y = np.exp(-x/1)*np.sin(50 * 2.0*np.pi*x)

    #Fourier
    yf = scipy.fftpack.fft(y)
    xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
    yp = 2.0/N * np.abs(yf[:N//2])

    label = "N="+str(N)
    ax.plot(xf, yp, label=label)

ax.legend()
plt.show()

The result gives different peak amplitudes for the different N-cases. The more sample points are taken (higher N), the lower the fft peak value... So, if you take N really high, then the Fourier transform amplitude decreases?
I would expect that the peaks have the same amplitude because $y(x)$ remains the same function in all cases. If the function does not change, the Fourier transform does also not change...
How to obtain the correct peak amplitude, such that the peak value is the same in the three cases?


texts: Does a closed loop system with a reference of its previous iteration always reach a point in which it adventually loops, even if that loop is stable?

texts:
For example you have a $F(\cdot)$ which starts off taking some input $a$, which describes some set of objects and their states, and creating a output $F(a)$, which can be renamed $F_1$, that also describes a set of objects and their states, which is then input into $F$, to get $F_2 = F(F_1)$, or $F(F(a))$.

If this is repeated for an extensive period is it guaranteed that that the values of the individual iterations will be loop? For example $F(F_{x+n})=F(F_{x+2n})$ etc and $F(F_{x+m})=F(F_{x+2m})$ where $x$ is the start of the cycle and $m$ and $n$ are the frequency of the cycle.


texts:
Project: Face Detection

Description: I want to detect and crop a face in an image. The image is captured through webcam and only one face per image.I used OpenCV face detector, but I was not satisfied with the cropping. So, I started using STASM (http://www.milbo.users.sonic.net/stasm/) face landmark detector to crop the image.STASM uses OpenCV face detector to find face in an image and STASM locates landmarks in faces. In bad light conditions, the cropped image from STASM is not good as it is not exactly detecting the face alone.

1) I want to know any better algorithm for face detection. My main aim is to crop the face from an image.

2) currently I am using STASM for cropping. In bad light conditions or when in an image, if the whole or complete face (forehead to chin) is not captured, STASM cropping is not reliable (The output will be only eye or lips). And in my application, if there is no proper output from the stasm or if the face is not cropped prpoerly then i should reject the images. How to do that? So I am planning to validate the face in an image by finding the Eyes. If I am right in my approach, how to detect the eyes from the cropped image?


texts: Autocorrelation function of a triangular wave

texts:
I am interested to find the analytical expression for the autocorrelation function of a signal that comprises triangular pulses. I have followed the derivation in "Statistical Theory of Communication" by Lee (1960, p.16)
However, if I replicate this in MATLAB, then I do not obtain the autocorrelation function as derived here. The general shape is the same, but it looks as follows:

with the autocorrelation function dipping below zero. In this case, the triangular pulses are spaced such that each triangle is half of the period, with the remainder of each period being zero.
Why is this different from Lee's Derivation, and how can I go about adapting his technique to find the analytical form that matches the empirical result?
The code I used in MATLAB was as follows:
function [x, acf1, acf2] = ACFTriangle(leftbase,rightbase,height,gap,repeat)

x = [zeros(1,gap) (1:leftbase)*(height/leftbase) (rightbase-1:-1:1)*(height/rightbase) 0];

x = repmat(x,[1,repeat]);

t = 0:length(x)-1;

acf = xcov(x,'unbiased');

n = length(acf);

s = (n-1)/2 + 1;

acf1 = acf(s:n)/var(x);

acf2 = autocorr(x,length(x)-1);

x = [t; x]';

end

The reason for my question is that I have a very large number of time series in which there may be a hidden triangular pulse wave. If I have an analytical expression for the autocorrelation sequence, then I can identify their characteristics from the shape of the autocorrelation function.
New edit
Notwithstanding the excellent answer below, how would one go about modifying Lee's mathematical derivation to allow for removal of the means? I can attempt this doing the following:
\begin{equation}
\varphi_{11} = \frac{1}{T_1}\int_0^{b-\tau} \left(\frac{E_m}{b}t - m\right)\left(\frac{E_m}{b}\left(t+\tau\right) - m\right)\,dt
\end{equation}
where $m$ is the mean of the series.  However, following this through gives:
\begin{equation}
\varphi_{11} = \frac{E_m^2}{6b^2T_1}\left(\tau^3-3b^2\tau + 2b^3\right) + m^2\left(b-\tau\right) - \frac{mE_m}{bT_1}\left(b-\tau\right)^2 - \frac{mE_m}{bT_1}\tau\left(b-\tau\right)
\end{equation}
which clearly gives $\varphi_{11}=0$ for $\tau=b$.  Where is the error in this working?


texts: Fourier Transform Computer Vision Textbook

texts:
I am required extremely fast to fill the gap and learn Fourier Transform with application in Computer Vision.

It was easy to find mathematical aspects of Fourier Transform, but I am more interested in application in Computer Vision for a newbie with a lot of examples.

If you familiar with a textbook or have any other source of information about Fourier Transform well explained with examples in the field of Computer Vision, please share it with us.


texts: How can we do finite-step signal sampling in practical?

texts: The signal has infinite length. We cannot sample it with a whole length. How can we truncate a infinite signal and sample it in finite step? And then how can we analyze it.


texts: Why $H(A)=H(C)$ where $C$ is $A$ with an additional parity bit?

texts:
Let $A=\{00,01,10,11\}$ with equal probabilities for each symbol, and $B=\{0, 1\}$ be a parity generator such that
$$
b=\begin{cases}
0, & \text{if} \,\, a=00 \quad \text{or} \quad a=11 \\
1, & \text{if} \,\, a=01 \quad \text{or} \quad a=10
\end{cases}
$$
Now assume we transmit $(a_i,bj)$ where $a_i$ is a symbol in $A$ and $b_j$ is the parity bit associated with it. I calculated the entropy for $A$ as follows:
$$
H(A)=4 [4 \log_2(4)]=2
$$
and to calculate the entropy of the new symbols, call it alphabet $C=\{000,011,101,110\}$, we need the probabilities:
$$
p_c(0)=\mathbb{P}(a=00 \, \text{and} \, b=0)=\mathbb{P}(a=00)\mathbb{P}(b=0 \mid a=00)=\mathbb{P}(a=00)=1/4
$$
Similarly, $p_c(1)=p_c(2)=p_c(3)=p_c(0)=1/4$, the entropy of $C$ is
$$
H(C)=4 [4 \log_2(4)]=2=H(A)
$$
How is this the case? We have more bits per symbols in $C$.


texts: Using large DFT for OFDM to avoid the guard interval

texts:
In OFDM, we usually insert the guard interval to avoid the inter-symbol interference at the expense of achieved data rate. that guard interval can be either zeros or CP guard interval, which it must be longer than the channel delay.
I am wondering, what's about using long symbol to avoid that guard interval ? for example instead of having number of subcarriers $N= 64$ or $N = 1024$, we can use $N = 5120$, what's the impact of using a high number of $N$ ??


texts: System Properties; Linear, Causal, Time-Invariant, Stable?

texts: I know the answers of the below question but I dont know why, except linear, could you explain the rest? Small tick: correct, small tick with dash incorrect answer.


texts: It is said that sparse code multiple access (SCMA) has the advantage of higher shaping and coding gain over LDS-CDMA due to the multidimensional constellation. What exactly is the shaping gain and how multidimensional constellation  increases the shaping gain in SCMA?


texts:
I am given 256 samples at a sampling rate of 48k and am asked to find the frequency of dominate tones and the amplitude of the largest tone. The data looks like this when plotted:

I assume I must use a FFT transform to do this. I've tried this in matlab:

Fs = 48000;
x = dataSet; % 256 samples data set
xdft = fft(x);
maxAmp = max(abs(xdft));
%Not sure how I can grab the tones here
freq = 0:Fs / length(x):Fs/2;

Not really sure how to proceed, Any ideas?


texts: Hamming window in frequency domain and back

texts:
I can't figure out the scaling and slight artifact in fourth plot.

N             = 32;
window        = hamming(N);

A             = fft(window, 256);

A_ifft        = (N/2) * ifft(A)

subplot(411)
plot(window);
ylabel('w(n)')
title('Hamming window')

magnitude     = abs(fftshift(A));
subplot(412)
plot(magnitude);
ylabel('W(j\omega)')

response      = 20*log10(magnitude);
subplot(413)
plot(response);
ylabel('W(j\omega) dB')

subplot(414)
plot(abs(fftshift(A_ifft)));
ylabel('IFFT(W(j\omega))')


texts: I can't implement plateau limit histogram equalization algorithm using MATLAB

texts:
i tried to implement this research paper using MATLAB. Unfortunately, my result is not same with the original result.

Here is my try step by step:

histogram_of_image = imhist(input_image);
a = uint8(2); %according to research paper
beta = uint8(2); %according to research paper
for i=1:256
  modified_histogram(i)  = (log( histogram_of_image(i)+double(a)))^double(beta);
end

for i=1:256
    if(modified_histogram(i)~=0)
        sum1 = sum1 + (modified_histogram(i));
        cnt = cnt + 1;
    end
end
tcl = sum1/cnt;
clipped_histogram = zeros(1,256);
for i=1:256
    if((modified_histogram(i)) >= tcl)
    clipped_histogram(i) = tcl;
    else
    clipped_histogram(i) = (modified_histogram(i));
   end
end

PDa = zeros(1,256);
for i=1:256
 PDa(i) = clipped_histogram(i) / (sum(clipped_histogram));
end
CDa = zeros(1,256);  %create CDa in formula
CDa(1) = PDa(1) ;
for i=2:256
     CDa(i) =  PDa(i) +  CDa(i-1);
end
value_after_enhancement =  zeros(1,256);
for i=1:256
     value_after_enhancement(i) = (255 * CDa(i)); %maximum gray level is 255 and minimum grey level is 0
end
b=uint8(0);
output_image=zeros(width,height);
for i=1:width %map new values
    for j=1:height
            b = input_image(i,j);
            output_image(i,j) =   value_after_enhancement(b);
    end
end

and if you would like to see, here is my all code


texts:
I've generated this spectrogram using wavelet convolution.
The chirp is generated like so:
framerate = 10000
DURATION = 10
N = DURATION*framerate
k = np.arange(N)
chirp = sig.chirp(k, 0.01, N, .49)

Each row (linearly representing a frequency) is generated like so:
$$row_f = |IFFT(FFT(chirp) * FFT(wavelet_f))|$$
where $f$ represents the frequency.
When I use matplotlib's built in spectrogram tool, I get this image:

L = 100  # length of one segment
overlap = L/2  # overlap between segments
plt.figure(figsize = (10, 5))
plt.specgram(sounddata, NFFT=L, Fs=2, noverlap=overlap, sides='onesided')

In this plot, the line is seemingly cleaner - not blurring over time. Why is the line blurring? I thought wavelet convolution was supposed to have higher resolution than STFT.


texts:
I'm having a really weird issue with generating waveforms from frequencies C1 to C8.

When I generate sine waves, everything sounds good; no problems whatsoever. However, when I generate a square wave... BAM weird aliasing issues on the higher frequencies (around C6 to C8). I have no idea what I'm doing wrong at this point.

I'm going to take a wild guess here but I'm thinking it might be because square waves have a lot of harmonics so it's not producing the higher harmonics correctly when I reach a higher frequency... thus the really high harmonics are getting aliased. Would like someone to confirm with me/guide me to a solution to this. Is it possible to try to implement a filter or is there a better solution ?

Thanks!

I'm setting my sample rate to 44.1 kHz and here are my functions for waveform generation:

float generateSin {
    sample = sinf(_phs);

    _phs += _phs_incr;
    _phs = wrapPhase(_phs);
}

float generateSqr {
    sample = sinf(_phs);
    if (sample > 0) sample = 1;
    else if (sample < 0) sample = -1;

    _phs += _phs_incr;
    _phs = wrapPhase(_phs);
}


texts: How to produce a high-pass filter from a low-pass one?

texts:
I already have a low-pass filter and want to remove the DC offset of a signal with a high-pass one.  I guess you can subtract a low-pass filtered signal from the original signal and the result is a high-pass filtered signal?

My low-pass filter:

x = input, raw signal
y = output, low-pass filtered signal
n = current sample
n - 1 = previous sample

y(n) = y(n - 1) + (x(n) - y(n - 1)) / alpha

High-pass filter:

x = input, raw signal
y = input, low-pass filtered signal
z = output, high-pass filtered signal
n = current sample

z(n) = x(n) - y(n)

Is this last equation correct?


texts: STFT with varying window lengths (like the continuous wavelet transform)

texts:
I was wondering if it is possible to use windows of varying lengths when making a spectrogram based on the short-time Fourier transform (STFT). That is, for higher frequencies I would use shorter windows so that I get higher temporal resolution at the cost of poor spectral resolution. I understand this method will result in tradeoffs between time and frequency resolution just like in the continuous wavelet transform, as well as higher computational complexity.
I wanted to ask if there are any other disadvantages other than the ones mentioned above when the STFT is employed in such way since I could not find any literatures about it.
Thank you very much.


texts: Load and Play Wave File (`.wav`) on MATLAB

texts:
I created a wav file by "audacity" app on my mac and imported in MATLAB. This is the simple code I wrote to read the file. But as it reads, it plays so much slower than its actual speed. Does anybody know how to solve this problem? Thanks.

d = audioread('HelloWorld.wav');
sound(d)


texts: Finding when a part of audio starts and ends in another audio

texts:
I have two audio files in which a sentence is read (like singing a song) by two different people. So they have different lengths. They are just vocal, no instrument in it.

A1: Audio File 1
A2: Audio File 2
Sample sentence :  "Lorem ipsum dolor sit amet, ..."

I know the time (in decimal seconds) every word starts and ends in A1. And I need to find automatically that what time every word starts and ends in A2. So, how to find sound of a word (in A1) in another audio (A2) that has different duration (of word) and different voice?


texts:
I am currently working on a problem on filtering. In fact, let's suppose a signal that is expressed as :
$r(t) = z(t) + x_i(t)$, where $z(t)$ is a smooth function and $x_i(t)$ is white Gaussian noise.

The question is to find the optimal filter for eliminating the Gaussian noise $x_i(t)$ if nothing is known about the smooth signal $z(t)$?".

Sure thing, Kalman Filter has been lately solicited faced to problems that ask for optimal filters. However,in this situation, that Kalman filter is the best candidate since we do not know anything about z.

Thus, I think they are better candidates for this problem.


texts: transfer function and 'causal' signal - evaluate transfer function or use z-transform of input?

texts:
From my studying difference equations and transfer functions, I understand that when a complex exponential input $x[n]=z_1^n$ is applied to an LTI system with transfer function $H(z)$, determining the output involves evaluating the transfer function at $z=z_1$, i.e. $H(z_1)$ which is a complex quantity that simply scales the input signal to produce the output.  The textbook I have uses a sinusoid input to show the effect (the transfer function causes a magnitude and phase change to the input).

The answers to this question I asked a while ago made the point of recognizing the difference between $a^n$ and $a^n \cdot u[n]$.  All of that was in the context of two sided signals and the bilateral z-transform.

In something I found online, in the context of causal signals and the unilateral z-transform, there was a question posed: "if $x[n]=(1/2)^n \cdot u[n] \to y[n]=\delta [n-2],$ what is the output for $x[n]=cos(\frac{\pi}{3}n)?$",  where the author determined the output by evaluating $H(z)$ determined from the first relationship at $z=e^{jn \pi /3}$.  Since this is a 'causal' context, isn't that cosine really $x[n]=cos(n\frac{\pi}{3}) \cdot u[n]$ and hence requires taking the z-transform of the input?

When can you just evaluate $H(z)$ at the value of interest, and when do you have to take the (unilateral) z-transform of the input?


texts: What are three different techniques to convolve a signal

texts: I have studied convolution in time domain and frequency domain. But doesn't know the 3rd method and technique


texts: Splitting an audio signal into multiple audio signals based on frequency range in Librosa

texts:
Using the Librosa package in Python, how may I separate an audio signal into multiple audio signals based on frequency range?

I have a file music.mp3. I used HPSS to separate the music's harmonic and percussive components, which I was able to do with librosa.effects.hpss(). This was done so I can use the percussive component for a more reliable beat detection using Librosa's beat tracker.

However, I'd like to further separate the harmonic component of music.mp3 into two more separate signals bass.mp3, midrange.mp3 where:

bass.mp3 contains the 20-250 Hz frequencies of the harmonic component,
midrange.mp3 contains the 251-4000 Hz frequencies of the harmonic component.

I plan on doing this because when constructing a chromagram out of the harmonic component, it only gives me the twelve pitch classes for that component, where it coalesces all frequencies of the same pitch class (e.g. C2, C3, C4) into one bin (C). By separating the harmonic component further, I am able to get more divisions.

But how do I do this with Librosa?


texts: How to calculate the energy levels of frequency intervals of an audio signal?

texts: Using the Fourier transformation I can translate a signal from the time-domain to the frequency-domain. If I divide the frequency axis into intervals, how can I calculate the audios signals energy for each of these intervals?


texts: How to get a real signal back to a non-complex form after processing it in frequency domain (using fft)

texts:
I have an audio signal that is composed of 2 channels. The signal is all real numbers.

I transform this signal into frequency domain using fft, multiply it with a filter frequency response, and then transform it back into time domain.

My problem is that the result is a complex signal.

How can I convert this complex signal back into a real time signal.

The filter is not a real filter.


texts: Goertzel algorithm produces incorrect phase

texts:
I've implemented Goertzel algorithm according to the Wikipedia page (https://en.wikipedia.org/wiki/Goertzel_algorithm) and another page (http://www.mstarlabs.com/dsp/goertzel/goertzel.html), which are consistent with each other. I am then testing it using a naive implementation for a single DFT frequency (https://en.wikipedia.org/wiki/Discrete_Fourier_transform). I am getting correct amplitude, but incorrect phase. Please help me find the mistake in my code, I've spent a couple of days on it now, and can't see the problem...

/**
 * Complex imaginary number i.
 */
template<typename T>
constexpr auto COMPLEX_i = std::complex<T> (
  static_cast<T> (0.0), static_cast<T> (1.0));

/**
 * Convert frequency in Hertz to Omega (radian per sample).
 */
template<typename T>
T
hertz_to_omega (T hertz, T sample_rate)
{
  return hertz / sample_rate * static_cast<T> (2.0) * M_PI;
}

/**
 * Simple implementation of the Goertzel algorithm
 * for detecting a given omega frequency from a range
 * of raw signal values.
 */
template<typename T, typename Range>
std::complex<T>
goertzel (const Range &range, T omega)
{
  const auto coeff = static_cast<T> (2.0) * std::cos (omega);

  T s1 = 0.0;
  T s2 = 0.0;
  for (const auto value : range)
  {
    const auto s0 = coeff * s1 - s2 + value;
    s2 = s1;
    s1 = s0;
  }

  return std::complex<T> (
    std::cos (omega) * s1 - s2,
    std::sin (omega) * s1) /
    static_cast<T> (range.size ());
}

BOOST_AUTO_TEST_CASE( test_goertzel )
{
  using T = double;

  const T FREQ = 543.2;
  const T SAMPLE_RATE = 5678.9;
  const T OMEGA = hertz_to_omega (FREQ, SAMPLE_RATE);
  const T START_RADIANS = M_PI / 3.0;
  const size_t LENGTH = 20000;

  std::complex<T> dft_sum (static_cast <T> (0.0), static_cast<T> (0.0));
  std::vector<T> data;
  {
    T radians = START_RADIANS;
    for (size_t sample_index = 0; sample_index < LENGTH; ++sample_index)
    {
      const T value = std::cos (radians);
      dft_sum += value * std::exp (- COMPLEX_i<T>
        * OMEGA * static_cast<T> (sample_index));
      data.push_back (value);
      radians += OMEGA;
    }
  }

  const auto result1 = dft_sum / static_cast<T> (LENGTH);
  const auto result2 = goertzel<T> (data, OMEGA);

  std::cout << "DFT: " << result1 << ", ABS: " << std::abs (result1) << std::endl;
  std::cout << "GOE: " << result2 << ", ABS: " << std::abs (result2) << std::endl;
}

Running the test above produces the following output:

Running 1 test case...
DFT: (0.250009,0.433004), ABS: 0.499997
GOE: (0.114572,0.486693), ABS: 0.499997

The magnitude is detected correctly (it should be 0.5), but the phases of the naive DFT implementation and the Goertzel algorithm output don't match. Please help me find the mistake in either my single-bin DFT or Goertzel implementation.

EDIT:

It seems that DFT calculation is detecting the phase correctly, as atan(0.433/0.25)/180 = 0.333 (and I specified the phase of pi / 3 for the generated signal on which I am testing), so its my implementation of Goertzel algorithm which is incorrect.


texts:
Laplacian of Gaussian formula for 2d case is
$$\operatorname{LoG}(x,y) = \frac{1}{\pi\sigma^4}\left(\frac{x^2+y^2}{2\sigma^2} - 1\right)e^{-\frac{x^2+y^2}{2\sigma^2}},$$ in scale-space related processing of digital images, to make the Laplacian of Gaussian operator invariant to scales, it is always said to normalize $LoG$ by multiplying $\sigma^2$, that is
$$\operatorname{LoG}_\text{normalized}(x,y) = \sigma^2\cdot \operatorname{LoG}(x,y) = \frac{1}{\pi\sigma^2}\left(\frac{x^2+y^2}{2\sigma^2} - 1\right)e^{-\frac{x^2+y^2}{2\sigma^2}}.$$ I wonder why multiply by $\sigma^2$ not $\sigma^4$ or anything else?


texts: Amplitude Response at greater than half the sampling frequency

texts:
I am hoping to clear up some confusion I have. In a lab I am taking, we analyzed the amplitude response of a simple system. We found that as we increased the input signal frequency to greater than half the sampling frequency, the output signal began to flat line. The picture below show the corresponding Amplitude Response.

I am aware of the Nyquist theorem, and I thought that  the reason the output signal flatlined at frequencies greater than 8 kHz was because the sampling frequency was 16 kHz.

But during a second experiment of generating a sine wave, I discovered the concept of folding about the Nyquist frequency. Now in this case, increasing the generated signal frequency to greater than the Nyquist frequency caused the output signal to not decrease in amplitude, but to become symmetric about the Nyquist frequency.

E.g.

Generating a sine wave of 14 kHz produced an output of a sine wave of 2 kHz.

How are these two concepts related? It seems in the first example, going above the Nyquist frequency caused the output signal to flatline. In the other case, generating a signal above the Nyquist frequency didnt cause the output to flatline, but caused the output signal to only have a different output frequency.


texts: filter out an expected but not regular signal

texts:
I am stuck with a case where I need to get rid of a portion of sinusoidal signal that suddenly appears in my main signal (sine wave).

my total signal will look like:

$
s(n) =
\begin{cases}
A\sin(\frac{2\pi}{8000}*n + \theta) + b\sin(\omega*n) & \text{if } 4500 \leq n \leq 5500 \\
A\sin(\frac{2\pi}{8000}*n + \theta) & \text{elsewhere}
\end{cases}
$

in the given example, I have $8000$ samples which represents one period of the main sine signal, but at a certain range I have the interference of another(weaker) sine signal, that occurs and disappears abruptly (for almost 1000 samples out of total 8000).

My previous approach:

The strange signal used to appear for a very short range (around 50 samples out of 8000), and also a less amplitude (experimentally), so I used to pass all my signal through a Least-squares estimator to come out with the closest $A\sin(\alpha)+B\sin(\beta)$. It used to give a satisfying results. But for the current conditions, least squares is really helpless in front of the big impact of this [noise].

How I am trying to solve

Now, I am thinking of correlation, I am trying to gather as more information as possible about the (noise) signal through experiments, once I can expect the exact shape (freq, length, ...). I am thinking of applying a correlation between $s(n)$ and the designed $noise$. [still in process, but I would love to hear guidance and advise].

Here I attach an illustrative image of the case


texts: Feature matching of images without corners

texts:
I am trying to stitch aerial images using feature-based registration. Most image pairs are matched well, such as this:

Others are not matched due to lack of common features:

As you can see, the images contain almost no corners (only fuzzy edges).

Is there a way to register such images? Please note that I need homography so template matching is not possible. The images are taken by flying drone wiggling in the air, hence the images have perspective distortions and feature matching is therefore essential.

Please note I have tried many point features (Harris, SIFT, SURF, MOPS) with similar (bad) results.

My ideas:

Use some line-based (Hough transform?) or area-based features (MSER?)
Use more images or video and then recover homography from dense data using Horn-Schunck or Lucas-Kanade

Here you can download the source images - it contains two pairs of images, the easier pair and hard-to-align pair:

Download original images


texts: Beginner level : Help in understanding from paper how the log-likelihood term has been obtained

texts:
Based on document : Practical Approaches to Principal Component Analysis in the
Presence of Missing Values

The document explains probabilistic approach to principal component analysis using Maximum A - Posteriori (MAP) estimation. The way MAP is used is When we compnesate for the lack of information due to limited observations with an a priori preference on the parameters based on prior knowledge. The prior knowledge is given by the prior's distribution. Then,

$$\hat{\theta}_{MAP} = \arg \max  p\left(D\vert\theta\right)p(\theta)$$

We maximize over the logarithm.

$p(\mathbf{x}_n)$ is normal with zero mean and unknown variance, and $\mathbf{y}_n$ is a linear function of $\mathbf{x}_n$ plus independent Gaussian noise where $n = 1...N$ represents the number of data vectors where, each vector $\mathbf{y}_n \in R^d$ and vector $\mathbf{x}_n \in R^q$. The missing/ latent variables are $\mathbf{x}_n$ and the complete data is formed from the observations $\mathbf{y}_n,\mathbf{x}_n$.

(please note that I have denoted $m_i$ as $\mu_i$)

Problem 1: In Section 6.4 Eq(39) is the cost function for MAP is given.

  \begin{align}
C_{\rm MAP}=&\frac{1}{v_y}\sum_{ij\in O}\left(y_{ij}-\mathbf W_i^{\rm T}\mathbf x_j-m_i\right)^2+N\log 2\pi v_y+\frac{1}{v_m}\sum_{i=1}^dm_i^2+d\log 2\pi v_m\\
&+\sum_{k=1}^c\left[\frac{1}{v_{w,k}}\sum_{i=1}^d w_{ik}^2 +d\log 2\pi v_{w,k}+\sum_{j=1}^n x_{kj}^2 +n\log 2\pi \right].\qquad (39)
\end{align}

Can somebody please help in making me understand how this cost function comes ? The priors for $X,m$ are from normal distribution (please refer Eq(19-20))
Problem 2: I am unable to follow how the last term in Eq(39) has been obtained. The first 4 terms are obtained from taking the logarith of these 2 terms : $\prod p(y\vert x) \prod p(m)$

Also, There is another expression $p(d_{ij}) = 1/N \exp(-{|| d_{ij} - dist_{ij}||}^2/\sigma^2_d) $ where $d_{ij}$ is the estimated distance between latent points $\mathbf{\hat{x}}_{i}$ and $\mathbf{\hat{x}}_{j}$ and $dist_{ij}$ is the true known distance between the corresponding high dimensional points $\mathbf{y}_{i}, \mathbf{y}_{j}$
$dist_{ij}$ can be expressed as :

$$dist_{ij} = (\mathbf{y}_{i} - \mathbf{y}_{j}){(\mathbf{y}_{i} - \mathbf{y}_{j})}^T = (\mathbf{W}^T \mathbf{x}_{i} - \mathbf{W}^T \mathbf{x}_{j} ){(\mathbf{W}^T \mathbf{x}_{i} - \mathbf{W}^T \mathbf{x}_{j} )}^T$$

$$ = {(\mathbf{x}_{i} - \mathbf{x}_{j} )}^T \mathbf{WW^T}(\mathbf{x}_{i} - \mathbf{x}_{j} )$$

So, $dist$ is a function of the latent variable. I was thinking if the prior $p(X)$ can be included with this but I don't know the rule for doing so.

How can I include this pdf as a prior or some other way into the likelihood expression?


texts: What are the frequency components of a sampled signal after low-pass filtering?

texts:
A continuous signal $x_a(t)$ is a linear combination of sinusoids of frequency 250 Hz, 450 Hz, 1 kHz, 2.75 kHz and 4.05 kHz. The signal $x_a(t)$ is sampled at $f_s=$1.5 kHz, and the resulting digital signal is converted into analog using a reconstruction filter  with a low pass filter with cutoff frequency of 750 Hz. The resulting signals is $y_a(t)$.

The question is, which frequency components will appear in $y_a(t)$'s spectrum? I tried writing a Matlab program where I simulated $x_a(t)$ and then filtered it with a Butterworth low pass filter. I then looked at its frequency spectrum and a component of 500 Hz appeared. I showed this to my professor but he said it was wrong, and he didn't explain why.

So I'm asking you to tell me what is the real answer and why.

Thanks a lot :)


texts: Removing marker pen artifacts from scanned microscope slide images

texts:
I am developing an image processing algorithm / plugin for Fiji (ImageJ) that processes scanned microscope slide images (20x magnification) and extracts features from the slides.

The slides contain histological tumor sections that are stained immunohistochemically with an antibody that identifies nuclei positive with the Ki-67 protein, along with a counter-stain that identifies all nuclei. The so called proliferation index, which measures how fast the cells are dividing, is then calculated as the ratio of positive-to-negative nuclei. For the feature extraction, I basically need to extract two color components from the RGB image; the one corresponding to the positive nuclei (diaminobenzidine; brown color) and the other corresponding the negative cells (methyl green; light green color).

I am using Fiji's Color Deconvolution plugin (http://www.mecourse.com/landinig/software/cdeconv/cdeconv.html) for separating the two stains. This works remarkably well, except when the slide has been drawn onto with a marker pen. Pathologists sometimes mark regions of interest on the slide with a marker pen so that it would be easier to inspect the slide later on. My data set contains slides that have been drawn onto with a red/blue/green or a black marker. The problem is that the presence of colors other than the two stains fails the color deconvolution algorithm totally.

 Above is an example of a slide at 20x magnification showing the edge of a line drawn with a black marker pen on the top-half of the image. The lower part of the image shows the brown positive nuclei and the light green negative nuclei.

Is there any chance for filtering out these marker pen artifacts? The effect is similar to looking through a welder's mask, and I'd wan't to somehow subtract the mask from the image. If I can't filter out the marker, I'll probably just need to manually crop out the portions having these artifacts.

Edit:

I've given high-pass filtering a try as suggested by Aaron. I used ImageJ's High-pass filter plugin, which creates a Gaussian blurred image from the source and then subtracts this blurred image from the original. The image is here (the radius / blur length or standard deviation of the Gaussian used was 20 pixels):

The artifact is almost gone and the nuclei are still recognizable, which is nice. The yellowish background color in the original (due to a calibration problem with the scanner) is also gone.

However, I'd also like to understand how this works. :) What does the Gaussian do to the RGB values of each pixel? How should I select the pixel radius for the Gaussian? Wikipedia's Gaussian Blur article wasn't all that clear to me. Are there other ways for implementing a high-pass filter for image data?


texts: feature extraction techniques for iris recognition

texts: I want to ask how I can divide feature extraction techniques to feature detectors and feature descriptors. I have big problem how to understand it. For example I can use Gabor filters (feature detector) to extract phase information (so the phase is feature) about image and then process this information to binary information like Daugman (I dont know if this is feature descriptor). I dont understand others approaches like Laplacion of Gaussian, SIFT, PCA, DCT and etc. Thanks for answer.


texts:
given:

$$ H(z) = \frac{4z(z-1)}{z-0.5} $$

I would say, when all poles are in the unit circles, the impulse response is right sided and causal.


texts: How to understand the energy centralizes at low frequency after making DCT on an image?What does energy mean?

texts: I've already know the low frequency represents a smooth area,and most part of an image doesn't have abrupt change which can explain natural image has more low frequency. But I just can't understand what does the term 'energy' mean.I would really appreciate if some can help me!


texts:
I was studying signal processing and I was frequently asked to verify if a certain covariance is possible for a given random function. I tried to check by verifying the property to prove it:
${γ}_{xx}(τ) = {γ}_{xx}(t+τ,t)$
but didn't quite got anything.
For exemple, if we suppose $τ_0$ a given constant, how would I prove that I can have the following covariance (or not)?
${γ}_{xx}(τ) = exp(-(τ- τ_0)^2)$


texts:
I am looking for a method to calculate the instantaneous pitch variation within a speech signal for time-warping application. Sample by sample frequency variation knowledge within a frame will give more accurate output in time-warping.

Till now, I read about a method based on Teager's energy operator. Please suggest me some other methods for the same.


texts: How do I attenuate a tidal signal in depth time series data

texts: I work with depth time series data from electronically tagged fish. When the fish spend time on the bottom we get a prominent tidal signal (approx 1.96 cycles per day). This interferes with our analysis of the data, so I have been experimenting with removing it using FFT - i.e. run FFT attenuate the bins around 1.96 cpd then inverse FFT. It nearly works, but the tidal signal varies over time (the time series can be up to a year long, with varying sample rates between 10 s and 2 mins). So in some parts the signal is removed, in others it is exaggerated. I have attempted to divide the time series into blocks and deal with each individually, but it doesn't seem to work. I work in C#. Any help and advice would be much appreciated!


texts: Removing known sequence from an audio sample

texts:
I have an audio recording of a conversation which is obscured by a movie playing on a television.

I have an audio recording of that same movie playing on that same television with no conversation in the background.  The recording device is identical as is its placement in relation to the speakers.  The volume is close to identical.

How do I remove the audio from the television so that I can unobscure the conversation?

Obviously ambient noise is different and system/thermal noise in the recording device will be different.  The phase of the two audio recordings will be different as there is some time error bounded by $1/F_s$.  ($F_s = 44100\textrm{ Hz}$)

I was unsuccessful in attempting to find a time offset between the two signals to perfectly align the audio.  I loaded the files in MATLAB and took the sum of the absolute difference between the two samples for multiple offsets and then selected the one with the lowest error.

The reconstructed signal being the difference between the shifted reference audio and the conversation.

 reconstructed_signal = conversation - reference_audio(offset:(length(conversation)-offset+1));

This time domain approach has not worked.  I am considering that I should do a difference of the frequency data of the two samples with a sliding window.

Do you think that this would be a good approach?


texts: Common methods of transmitting digital data through an extremely noisy audio signal?

texts:
I am developing software for transmitting data from a computer to a mobile device using audio. Specifically, from the computer's speaker to the device's microphone. (I am developing both the sender and the receiver software.) This is one channel, and the sample rate is a standard 44.1Khz. It needs to be VERY resilient to noise - it will be used in environments with street noise, people talking, background music, etc., which may be 10dB louder than the signal itself.

Knowing almost nothing about DSP, thus far I have apparently implemented frequency-division multiplexed on-off keying, without realizing it until now :) Basically each of 96 frequencies corresponds to a bit, and on the receiving end, if the frequency is present, the bit is on. I then use some custom error correction code to whittle those 96 bits down to 8 very reliable bits.

My system works very well and is very robust, but unfortunately has an EXTREMELY low data rate (~80 bits/sec). I have some ideas about how to improve bandwidth, like using more than 96 frequencies, and using the phase component to encode additional data. Also it appears I should look into frequency-shift keying and phase-shift keying.

So - my question is simply this:

What are the options for encoding data in a noisy broadcast audio signal, and what kind of bandwidth and error rates can be expected with each technique?


texts:
Given a filter used to shape the digital signal, $p(x)$, and given that we do not want the filter combination to cause any ISI, what "matched" filter, $q(x)$ will maximise the SNR?

Matched filters are used in digital communications to maximize the signal to noise ratio.  Often a root-raised-cosine filter is used to shape the signal, since it is bounded in frequency space and the same filter can be applied to the received signal to improve the signal-to-noise ratio (SNR) without causing inter-symbol-interference (ISI).

However if a less optimum filter is used to shape the signal, then using the same filter at the receiver can introduce ISI.  It is not immediately obvious what the best choice of filter at the receiving end is.

My understanding is that the SNR is maximized by maximizing $\int{p(x)q(x)dx}$, so I want to maximize this while satisfying the constraint that filters cause no ISI ($p(x)*q(x) = 0$ for $x=kT$, $k$ is an integer, $T$ is symbol width).

Presumably one could do this by solving a Euler-Lagrange equation with some lagrange multipliers for the constraints.  Is there an easier way, or am I making a mistake, or going in the wrong direction?


texts: Question about the Haar Wavelet transform

texts:
I'm trying to detect the presence of a sinusoid in 1/f noise conditions using the STFT and the DWT with Haar wavelets. I find this interesting phenomenon that I'm not able to explain (see plot at 1)

My sampling rate is 1200 Hz and the graphs show the one sided response (0-600 Hz). I am using 2 scales for the Wavelet Transform.

My question is, why in the FFT method can I clearly see the blips in about the correct frequency as my simulated data (180 Hz), while in the WT method, the increase in energy is spread visible across all three bands. Shouldn't I see a bright line only in the highpass segment of the second scale?


texts: Is this system an FIR filter or an IIR filter?

texts:
Consider this system:

 y[n]=y[n-1] + x[n]

It can be verified easily that this system has a finite impulse response (by putting x[n]= delta[n]  i.e impulse, and evaluating impulse response h[n]=y[n]). Yet for a unit step input U[n], this system is an unstable system.  On wikipedia it is mentioned that FIR filters are always stable

Now my question is , is this an FIR filter or an IIR filter?


texts: How can I load a wav file, split it in frames, apply different filters do different frames and then recombine the filtered frames into a single file?

texts:
I'm working on a wav audio file using MATLAB. What i would like to do is:

load a wav file and divide it into frames of the same dimension
filter these frames (with a different filter for each different frame) using, for instance, iirnotch
put all the filtered frames together again without having discontinuities between them

I already wrote a script using the overlap-and-add technique using a "Hann" window and it works quite well but between the different frames I'm not able to avoid the discontinuities so that when I listen to the wav file (made by adding these frames together again) I hear an annoying glitch disturbance.

Anyone know how to avoid this problem?


texts: What is the Fastest Way Ever to Do 2D Correlation?

texts:
Here is what i gathered so far:

Any sliding window classification, image filtering or similar can be fastly done by a FFT (flip the signal and do convolution).
Also, if the template/filter kernel is separable it is possible to do fast correlation by simply separating into multiple kernels and applying then sequentialy.
I think it is possible to do fast correlation by using integral images, but I am not sure how.

The question is: is FFT the fastest way? If now, which is? Consider both GPU and CPU implementations

If it is not: is there algorithm capable of separating a non-separable kernel (in aproximated way, logically)?

Context:

My filter is an Inner Product Detector, which is found by

$$
h = R^{-1}\mu_1
$$
where $ R^{-1} $ is the autocorrelation matrix and $\mu_1$ is the mean of a given class.

The matrix $ R $ is hermitian. And can be calculated by

$$
\sum x^Tx
$$

where $x$ are d-dimensional vectors and $[.]^T$ denotes the transpose.

Also: kernels are squares of size 21,23,25,27 or 35
Images are 128x128


texts: What is Causing a Noise Signal in USB Audio Recording?

texts:
I am experiencing the following peculiar noise signal when recording from my Behringer Xenyx Control 2 over USB.

There is frequency at around 2293 hz as show:

And over time, it produces the following plot:

What on earth is doing this?


texts: How to solve this math problem regarding Histogram matching?

texts:
Please help solve the problem in the picture. I am confused due to the piecewise nature of the density functions. If it was a simple curved PDF, I would simply integrate over full range and use the single CDF expressions. But, for this piecewise case, in each section I will have different expressions.
P.S.--This is NOT a homework question. It is from a past year exam.


texts: How to increase the spectral resolution?

texts:
My question is about the spectral resolution of a discrete signal. Each sample of my signal is made up with 2^n frames sampled at 44.1 kHz.

So, when I want to know the spectral resolution, I calculate : 44100/number_of_frames. With 2048 frames, my spectral resolution is around 20Hz. But, when I take a look to the bands filtered by an equalizer, the band with in the low frequencies are around only 5 Hz (16 Hz -> 20 Hz -> 25 Hz -> 32 Hz...). How is it possible ?

I thought about using "zero padding", but, even if it will help me to have a better location of each peak of the the spectrum analyze, this method don't magically increase the spectral resolution.

I also thought about increasing the number of frames analyzed. But, to get a spectral resolution of 5 Hz with a signal sampled at 44.1 kHz, I would need 8192 frames and it would represent 185 ms. It's very far from a pseudo real-time analysis and a singer who would listen his voice after this analyze while he is singing would hear this "delay".

So, what is the solution ?

Thank you for all your reply


texts: Proving that a product of matrices invertible

texts:
Given $R_x$ a Positive Definite (PD) covariance matrix of size $M\times M$ and $C$ a full rank $M \times N$ matrix, I want to prove that $C^* R_x^{-1} C$ is invertible to derive the Linearly Constrained Minimum Variance Beamforming.

My ideas so far:

Since $R_x$ commutes with its adjoint, it can be written using the eigendecomposition $R_x =  U \Lambda U^*$
$R_x$ is PD, $R_x^{-1}$ is also PD since $\lambda_i > 0 \Longrightarrow \lambda_i^{-1} > 0$
Matrix $U$ defines an orthonormal basis and $C$ is full rank. Can view $U^* C$ as a projection onto this basis.


texts: Best method to extract phase shift between 2 sinosoids, from data provided

texts:
I have been asking around if the way I was extracting phase shift (lag) was correct, and I ran into some trouble.

So in general, given 2 arrays of data of the same length, representing:

the input sinusoid and
the response, also a sinusoid.

So not knowing the functions that produced the input (but knowing the frequency, starting phase of the input, fs, and so on), what is the best way to find the phase lag?

should I do it "manually" by looking at the graphs and finding the phase lag at each component?
should I use FFT, like I have been trying?

any suggestions, links, books are welcome.


texts: Determine Parameters of gaussian distribution to fit a certain Image-Filter

texts:
I'd like to know what the parameters for the two-dimensional gaussian distributions are, which create this (I think) quite common image-processing filter:
$$
h_\sigma(x,y)=
\left[
\matrix
{
1&2&1\\
2&4&2\\
1&2&1
}
\right]
$$
The two-dimensional gaussian equation looks like:
$h_\sigma(x,y)= A\cdot e^{-\frac{x^2+y^2}{2\sigma^2}}$, where $A = \frac{1}{{2\pi}\sigma^2}$, but I think $A$ is not important as it is only a factor, isn't it?

With $\sigma=1$ I get:
$$h_{\sigma=1}(x,y)=
\left[
\matrix
{
h(-1,-1)=0.3679&h(0,-1)=0.6065&h(1,-1)=0.3679\\
h(-1,0)=0.6065&h(0,0)=1.0&h(1,0)=0.6065\\
h(-1,1)=0.3679&h(0,1)=0.6065&h(1,1)=0.3679
}
\right]
$$
As you see, I think the entries in the filter-matrix correspond to sampling values of the gaussian function from [-1|-1] to [1|1]. I hope that is at least right.

Multiplied $h_{\sigma=1}(x,y)$ with $4$, and then rounded gives me my desired filter, but I don't think that is the correct way, is it?

$$\mbox{round}(4\cdot h_{\sigma=1}(x,y))=
\mbox{round}(4\cdot
\left[
\matrix
{
0.3679&0.6065&0.3679\\
0.6065&1.0&0.6065\\
0.3679&0.6065&0.3679
}
\right]
)=
$$
$$
h_{\sigma=1}(x,y)=
\mbox{round}(
\left[
\matrix
{
1.4716&2.4260&1.4716\\
2.4260&4.0&2.4260\\
1.4716&2.4260&1.4716
}
\right]
)=
$$
$$
\left[
\matrix
{
1&2&1\\
2&4&2\\
1&2&1
}
\right]
=h_\sigma(x,y)
$$

There must be a better solution, or mathematical explanation for creating that matrix? I especially dislike the rounding operation, as the values are all close to $.5$ and not to $.0$ which should results in a huge error...

My question in other words: How do I determine the filter from scratch...


texts: where is the matched filter being used exactly in Rx PHY?

texts:
Recently I learned about the matched filter and its usefulness as it maximizes SNR at time t0.
So I looked up online to see how this amazing filter is being used exactly in 802.11a Rx PHY, but couldn't really point out which blocks in Rx PHY the matched filters are used exactly.

I suspect that matched filter can be used as a demodulator, but if one deploys 64-QAM or 128-QAM, do we still use matched filter for the demodulator in this case as wel??

so I have two questions :
Q1: Where in Rx PHY is the matched filter being used exactly?
Q2: If matched filters are used in demodulator as I expect, how is it implemented for a complex demodulation scheme (i.e., 64-QAM).


texts: Autocorrelation of filtered cyclostationary random process

texts:
Background:

When a wide-sense stationary (WSS) random process $x(t)$ is passed through a linear, time-invariant (LTI) filter with impulse response $h(t)$ to produce $y(t)$, the following relationship holds:
$$
S_y(f) = |H(f)|^2 S_x(f) = H^*(f) H(f) S_x(f),
$$
where $S_x, S_y$ are the power spectral densitities (PSD) of $x$, and $y$, respectively, and $H(f)$ is the Fourier transform of $h(t)$ (ordinary frequency, unitary).  In the time-domain, this relationship can be shown to be:
$$
R_y(\tau) = h^*(-\tau) * h(\tau) * R_x(\tau),
$$
where $R_x,R_y$ are the autocorrelations of $x$, and $y$, respectively.

Question:

If $v(t)$ is a wide-sense cyclostationary process with autocorrelation $R_v(t, \tau)$, is cyclic with period $T_0$, and is passed through the same linear filter to produce $w(t)$:
$$
w(t) = h(t) * v(t),
$$
will $w(t)$ also be cyclostationary with cyclic period $T_0$? (Assuming $h(t)$ is BIBO stable.) If so, is there a "simple" relationship that expresses $R_w(t, \tau)$ in terms of $h$ and $R_v$?

EDIT:

For a wide-sense cyclostationary process $v(t)$, the cyclic autocorrelation $R_v^{n/T_0} (\tau)$ is defined as:
$$
R_v^{n/T_0} (\tau) = \int_{-T_0/2}^{T_0/2} R_v(t, \tau) e^{-j n 2\pi t/T_0} dt
$$

A relationship between $R_v^{n/T_0} (\tau)$ and $R_w^{n/T_0} $ would also be acceptable for this question.


texts: Using DFT Rotation for Frequency Down-Mixing

texts: Is it possible to perform Frequency Down-Mixing (sometimes referred to as Down-Conversion) by simply rotating the DFT sequence? If so what are the Advantages and Disadvantages of this method?


texts: Method for calculating musical note-based magnitude spectrum

texts:
My specific question is: is there a standard way to do the kind of analysis I describe below? If not, is there a useful method that I have missed?

I'm exploring ways to produce a spectral histogram of audio signal magnitude based on musical notes. I wish to calculate a spectrum spanning 10 octaves, or 120 notes. This covers a frequency range from ~20Hz to ~21,000Hz. The sample rate will be at least 44,100Hz.

The final desired result is a numerical value per note, that I can then plot on an axis of note bins. This is essentially a logarithmic x axis, with non-linear note bins. I wish to do this in realtime, on a desktop computer, and display visually at ~10Hz update rate or better.

I initially looked at the DFT (via FFT), however this has some drawbacks - the size of the DFT needs to be large enough to provide sufficient resolution at the lowest frequencies. E.g. to get ~1Hz resolution at 44100 Hz I will need to use an FFT of size 65536, which is longer than 1/10th of a second of signal (for 10Hz update). Secondly, the frequency bins for the DFT are linearly spaced in frequency, so there needs to be a warped mapping from these linear frequency bins to the non-linear note bins. In some cases many frequency bins will map to a single note bin. I understand that the shape of the frequency bins is related to the window function, so I imagine I would use something like sinc interpolation to collect frequency bins into corresponding note bins.

I've also come across some other options:

Someone claimed to have found a way to tweak the DFT to use note frequencies as the basis functions, however I wasn't able to understand how they did this and there was no math provided to back it up. Is this possible? Can the set of basis functions for the DFT be changed to match note frequencies, and if so what are the implications for scaling, phase (not so important in my application). I'm a layman but my feeling is that if the basis functions are not orthogonal then there's going to be overlap which could result in some terms counting energy more than once.

Another suggestion involves applying the Goertzel algorithm 120 times. However based on the rule-of-thumb on Wikipedia, this would be significantly more expensive than using the FFT above.

Another idea involved applying a much smaller FFT to a single octave, with filtering and downsampling to shift each octave down. This would result in 12 smaller FFTs and 11 frequency-shift/filter operations. This may be more efficient than the large FFT above.

One final point - the system can be fine-tuned in realtime, so "A440" may not actually be exactly 440Hz and may shift slightly by up to half a semitone in either direction. Therefore I need a method that can be "tuned" accordingly, or is precise enough to correctly capture energy within +/- half-semitone of a note.


texts: Why does a strange frequency appear in this image?

texts:
Here is an image of a filter bank spectrograph, where the main signal is at 2.5kHz and the high signal is at 21.5khz,

What is the high frequency signal? I doubt that birds can make frequencies above 18kHz, so it is probably an effect from the recording equipment? Is it an issue from generic FIR BPF at high frequencies?


texts: Discrete-time model of 2-user CDMA system, recovered information at receiver using polarity detection?

texts:
I have discrete-time model of 2-user CDMA system with: $b_1=-1, b_2=+1, \alpha_1, \alpha_2, s_1, s_2$ and $n$ white noise.

What are the recovered infromation bits for user 1 and user 2 respectively if polarity detection is used?
If the optimum MLE is used, what are the recovered infromation bits?

I tried to solve but stucked with polarity.


texts:
I have the following problem with convolution operators. I know that if I have a filter $g$, then the adjoint filter would be given by $h(x,y)=g(-x,-y)$, so that $\langle g * u,v \rangle = \langle  u,f *v\rangle $ holds for the $\mathbb{R}^2$ scalar product.

I also know that $Ff(-\omega) = \overline{Ff(\omega)}$ for real $f$. So it seems I should hold that

$ h(x,y) = F^{-1}\overline{Fg} = F^{-1}Fg(-x,-y) = g(-x,-y)$

However when I compare

g = ifft2(fft2(g))

( for some gaussian kernel, e.g matlab: g = padarray(fspecial('gaussian',11,2),[5,5]) )

and

h = ifft2(conj(fft2(g))

I get the expected result for g

But when I compute h, which should be equal to g because g is gaussian

the result is shifted by (1,1)

which is somewhat unexpected.
I wonder if there is some subtlety in the discrete <-> continuous conversion that I have missed ?


texts: Spherical equivalent of Nyquist frequency

texts:
Let $\phi$ be a scalar function defined on the surface of a sphere. I have samples of $\phi$ at various locations on the sphere. I want to apply a spherical harmonic transform. I know that $\phi$ is 'band-limited' in the sense that it can be accurately represented by a truncated spherical harmonic expansion, which only includes terms up to a certain angular order, $l_{max}$. (Using this notation: http://mathworld.wolfram.com/SphericalHarmonic.html.)

The sampling is much denser than I require to accurately perform the spherical harmonic transform. I would like to discard some of the samples before doing the transform, so that the computation will be faster. How many samples do I need (in terms of $l_{max}$) to ensure 'good' results?

I realise that the answer will depend on the definition of 'good' and the distribution of the sample points (in my case, the samples are randomly distributed). However I am looking for a result like Nyquist's theorem from which I can make a safe estimate of how many samples to keep.


texts: Calculate the magnitude and phase of a signal at a particular frequency in python

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

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


texts:
I saw a video where this guy used a program to do a frequency analysis on a voice signal and a sawtooth wave (I'm assuming this was FFT). Then he saved the plots as images and combined them pixel by pixel through multiplication using photoshop. He then put this picture back into the program and it did the inverse transform, turning it back into a sound. He said that this was an implementation of a vocoder but the low quality made it hard for me to tell if he was correct.

If so, then couldn't one simply implement a vocoder as a convolution operation? You would just have your two signals then select a window from each of these of the same width at the same position, then perform the convolution operation on these two windows (and probably use a window function as well, such as Hanning). You would, of course, have to do this for every sample, so you would be doing this as many times as you have samples in your tracks (and the windows would sometimes lie partially outside of a track, so they would have to be zero padded).

This seems like it might work because the convolution theorem says convolution in the time domain is itemwise multiplication in the frequency domain, so if it doesn't implement a vocoder, it at least implements exactly what the man in the video was doing (at higher quality). And, I'm not just asking this question blindly, I actually tried it. I get a very cool voice effect but I'm not sure it's the same as a vocoder. In fact, it sounds nothing like a vocoder. If some kind individual could tell me exactly what is happening here, I would be very grateful.


texts: How increasing the Signal to Quantization noise increases the resolution of ADC

texts:
Inside ADC we have a quantizer, when the sampled signal is passed through it, the signal values get discretized. The number of discrete levels in which the signal is discretized is called the resolution of the quantizer. I read that if we oversample the input signal and then pass it through the quantizer and pass the output again through a low pass filter called decimation filter the signal to quantization noise (SQNR) of the quantizer improves resulting in better resolution of the quantizer.
How is this possible? True that the SQNR improves but how does that  affect the resolution of the Quantizer? Isn't the number of levels inside the quantizer fixed by its hardware and thus resolution of the quantizer i.e., the number of dicrete levels in which it breaks the input fixed?


texts: Introduction to fast algorithms for time–frequency analysis

texts: I had preliminary knowledge of digital signal processing from Oppenheim's Discrete-Time Signal Processing and is studying time-frequency analysis now. May someone suggest introductory reference (textbook, website, review paper...) for the fast algorithm aspects of time–frequency analysis?


texts: How can I trace a path around an image consisting of only black and white lines?

texts:
I'm looking for a way to trace a closed curve (i.e. a long list of pixel values all connected by lines in between) that approximates a bitmap image of black and white lines.

The image starts out as an arbitrary bitmap image. I've been using this one for testing.

Then it's converted to greyscale, passed through an edge detection filter, and finally converted to a true monochrome black/white bitmap. The end result of this preprocessing looks like this.

My goal is to trace a closed curve around the white pixels.

My naive approach was to make a list of all the white pixels in the image, then use a greedy travelling-salesman solver to draw a path around them. Which results in a path looking like this.

This works, there's a lot of artifacts but it works. But it's very slow. This relatively simple image takes a couple minutes to process. Any images more complex than this one need to be reduced in quality by a substantial amount before the time taken is reasonable. The resulting path is then almost unrecognizable.

When looking up algorithms for tracing images, I find a lot of things either built into image editing software (I need an algorithm that I can implement in Python), or overkill bitmap-to-vector algorithms designed to preserve color and line thickness.

All I need is a list of pixel values that trace out the lines.

There's got to be a better way than straight up solving travelling-salesman, right?


texts:
I am highly confused about the stability of $\cos(x[n])$.

If we provide a bounded input such as $x[n]=u[n]$, the output is bounded.

Now if we provide an unbounded signal $x[n]=\delta[n]$, the output oscillates between $-1$ and $1$. So I think it's marginally stable in this case.

Am I correct? I have this confusion because in my book it says it's unstable when $x[n]=\delta[n]$. So my questions are: Is the function BIBO stable and also is it stable for an unbounded input ?

Note: $u[n]$ and $\delta[n]$ are unit step and dirac delta functions respectively.

Correction: As mentioned in the comment, $x[n]=\delta[n]$ is bounded too.

But this points more strongly to $\cos(x[n])$ being stable, isn't it?


texts: USRP + GNU Radio as 1 GHz bandwidth spectrum analyzer

texts:
I have maybe not trivial but some interesting question.
I need to create 802.11ac channels spectrum analyzer with USRP.

Hardware: USRP B210
Software: GNU Radio 3.7.9
PC: Intel Core i3 i3210m 16 GB RAM.

According to documentation from: https://www.ettus.com/content/files/b200-b210_spec_sheet.pdf USRP has 56 MHz bandwidth, so FFT analysis greater than 56 MHz did it (by my opinion), by some multiple snapshots with 56 MHz.

During to experiments FFT with 1 GHz bandwidth refresh 1 per 3 seconds, so it can not be used for real time monitoring in that condition.

Of course it takes some processing time, but can somebody tell, it is possible to speedup FFT analysis on big bandwidth: 500 MHz, 1 GHz, at least to 1 refresh per 1 second, or it is not possible due to processing limits.

As my opinion, maybe need to investigate verilog code, that uploaded to USRP FPGA before first start, maybe it can be increased reference frequency with frequency multiplier, or not?

Marcus Müller, thanks for reply. I don't know how it can be enough detailed architecture explanation, but try to explain. user34533 = WhiteBear.
Sorry i use stack exchange first, and didn`t setup correct login from stack overflow. In future i will put posts from my regular account.

So:
USRP connected to the laptop via USB cable. In laptop installed Linux Xubuntu and GNU Radio Companion software. GRC interact with USRP via USRP source abstraction block, that powered by official UHD (USRP hardware driver), provided by Ettus Research also as open software. FFT block - standard GNU Radio block, and did by laptop. Bandwidth that displayed on FFT directly related to sample rate of USRP, so maybe i didn`t correct inform You. So when it setup to large values (greater that 400 MHz) FFT became continuously slow and slow refresh the spectrum.

Flowcharts and settings of blocks You can find on the images, ready to provide additional information. Thanks a lot for reply, hope, You can help


texts: The computation procedure problem of calculating the BER

texts:
Here is the processing of calculating the BER average

\begin{align}
P_e &=\int^\infty_0Q(\sqrt{2uSNR})e^{-u}du \\
 & =\frac{1}{\sqrt{2\pi}}\int^{\infty}_0\int^\infty_{\sqrt{2uSNR}}e^{-\frac{t^2}{2}}e^{-u}dtdu \\
& =\frac{1}{\sqrt{2\pi}}\int^{\infty}_0\int^\frac{t^2}{2SNR}_{0}e^{-\frac{t^2}{2}}e^{-u}dtdu \\
& =\frac{1}{\sqrt{2\pi}}\int^{\infty}_0(1-e^{-\frac{t^2}{2SNR}})e^{-\frac{t^2}{2}}dt \\
& =\frac{1}{2}(1-\sqrt{\frac{SNR}{1+SNR}})
\end{align}

channel state $h$ is a complex gaussian r.b.$CN(0,1)$,$SNR=E_s/N_0$,the BER of received symbols is $P_{e|h}=Q(\sqrt{2|h|^2SNR})$

I don't really understand these computation procedure step by step,for example.

Why is $\int^\infty_{\sqrt{2uSNR}}$ become $\int^\frac{t^2}{2SNR}_{0}$ ?
Why is $\int^\frac{t^2}{2SNR}_{0}e^{-u}du$ become $(1-e^{-\frac{t^2}{2SNR}})$ ?

Why is $\frac{1}{\sqrt{2\pi}}\int^{\infty}_0(1-e^{-\frac{t^2}{2SNR}})e^{-\frac{t^2}{2}}dt=\frac{1}{2}(1-\sqrt{\frac{SNR}{1+SNR}})$ ?

Can anyone explain them and teach me?