texts:
Pretty much a newbie here, but would like to understand a pretty basic concept.

We all know that a quick and dirty way of calculating signal power of a part of a discrete sampled system (such as audio signals) is the age-old power formula for RMS power, where we simply break the signal into chunks of N samples and sum the squares of the amplitudes, then divide by the number of samples

$$
P_\mathrm{x} = \dfrac{1}{N}\sum_{n=1}^{N}|x[n]|^2
$$

So let's say I have 10ms samples of two completely uncorrelated signals at 48khz (480 samples), and I have calculated the RMS power of these samples using the above formula and have come up with P1 and P2, the power of my small segment of each signal.

Then if I sum the two signals (simple software mixer function), what is the resultant power of the two summed signals?  Is it simply P1 + P2?

If so, do I simply keep going?  So if I had four signals with P1, P2, P3 and P4, if I sum all four of them together, is the resultant RMS power P1+P2+P3+P4?


texts: Teager-Kaiser Operator vs. Hilbert Transform

texts:
Since a couple of months I started working on the extraction (estimation) of signal frequency and amplitude components by means of two different time-frequency approaches, namely the Hilbert transform and the Teager-Kaiser energy operator.

I tested both methods on standard signals, such as chirps, sine, cosine, etc. This preliminary analysis seemed to highlight Teager-Kaiser has better resolution in both frequency and amplitude estimation.

Afterwards I applied the two methods on some acceleration signal derived from a tool simulating the dynamic behavior of a given system. Surprisingly enough, Hilbert transform provided more reliable results in terms of frequency estimation. Teager-Kaiser operator shows some high-frequency content estimation, which is hardly expectable.

I developed the two techniques in Simulink, as visible here:

The results on a synthetic acceleration signal are here:

What did I miss in the design of the Teager-Kaiser operator? Is it any possible to avoid those "false" frequency estimates (those with the arrows)?


texts: Determining whether or not a system is BIBO stable

texts:
I'm trying to determine whether or not a system with impulse response

$$
h(t) = \sum_{n=-\infty}^{\infty} \delta(t-2n)
$$

is BIBO stable. I haven't touched this material for a very long time -- could anyone lend a helping hand? I recall needing to show that

$$
\int_{-\infty}^{\infty} |h(t)| \, dt < \infty,
$$

but $h(t)$ is defined in a very weird way, and calculating the integral is proving to be difficult.

EDIT: Does the following argument seem reasonable? Since were summing an infinite number of $\delta$ functions, each will have an argument of $0$ for some unique value of $t$. Thus, there are infinitely many values of $t$ for which $\delta(t-2n)=1$, so the integral diverges.


texts: Book suggestion for upsampling and downsampling

texts:
I want to ask for a  good textbook or online material  which explains
in detail:

-the time and frequency domain mathematics
behind  upsampling and downsampling.
-easy to follow

Note: Oppenheim's discrete-time signal processing book is hard to follow.


texts: Different "zoom" after reconstruction from Radon transform using Fourier Slice Theorem

texts:
I have implemented the procedure to reconstruct an image from its Radon transform involving the Fourier Slice Theorem (FST). After application of the FST, the data is transformed from polar coordinates to cartesian, by regridding.
It works fine, except the reconstructed image is "zoomed out" (see picture). What causes this difference, and how can I resolve it, so that the phantom in the original, an reconstruction are exactly the same size?

Code:
N = 128;
f = phantom(N);

% Radon Transform
theta = 0:1:180;
[g, xp] = radon(f, theta);

% Fourier Slice Theorem
ghatpol = fftshift(fft(ifftshift(g)));

% Transform to cartesian coordinates
phi = -theta*2*pi/360;
[pX, pY] = pol2cart(phi, xp);
x = -N/2:1:N/2-1;
[X, Y] = meshgrid(x, x);
fhatcart = griddata(pX, pY, ghatpol, X, Y);

frec = ifftshift(ifft2(fhatcart));

figure,
subplot(121), imagesc(f)
subplot(122), imagesc(abs(frec))


texts: Fourier Transform gives unexpected results: signal reversal and conjugation

texts:
As far as I understand the math you can reverse a real valued signal by Fourier transformation, taking the complex conjugate of the result and inverse Fourier transformation, i.e.
\begin{equation}
\tag{1}
f(-x) = \mathcal{F}^{-1}(\overline{\mathcal{F}(f(x))})
\end{equation}
where $\mathcal{F}$ is the fourier transformation operation and the overline denotes complex conjugation. This all makes sense and from a mathematical perspective this isn't very useful. However, from a DSP perspective its very useful because the autocorrelation of a signal can be conveniently cast as the convolution of a signal with a time-reversed version of itself. If you're using an FFT convolution it means that you only have to calculate one FFT of the signal.

When I tried to confirm that equation (1) was working in python I found that the result is shifted by one sample. Here's the example code:

import numpy as np
from numpy.fft import rfft, irfft
x = np.arange(128)
y = np.sin(2 * np.pi * x / 128 * 4 + np.pi / 3)
y_rev = y[::-1]
y_rev_fft = irfft(conj(rfft(y)))
plt.plot(y, label="Original Signal")
plt.plot(y_rev, label="Straight Reverse")
plt.plot(y_rev_fft, label="FFT Reverse")

Is this the expected result? If so, can the conjugate still be used in the FFT correlation algorithm I described?

Also, is this question better suited to another stackexchange site?


texts:
I'm working with AEC of Speex. The algorithm is based on the MDF adaptive filter + an adaptive learning rate. I'm using it like a ANC and it works very well.
Does anybody have some material, as block scheme, data flow diagram of the AEC of Speex.
I have read the documentation but it is not very useful for the dsp concepts.
thank you so much
m


texts: Why do the two methods give different answers for the Fourier transform of $Y = \cos(\omega_0 t + \phi)$?

texts: Why do the following two methods give different answers (or are they the same) for the Fourier transform of $Y = \cos(\omega_0 t + \phi)$, with respect to $t \to \omega$ ?


texts: Time-variance and causality of $y(t)=\sum_{k=-\infty}^{+\infty}x(t)\delta(t-kT)$

texts:
Is the following system time invariant and/or causal?

  $$y(t)=\sum_{k=-\infty}^{+\infty}x(t)\delta(t-kT)$$

I think that it isn't time invariant and it is causal, and my thought is based on the fact that this is actually a sampling system that samples on $T-$intervals the input function. Is this correct and if yes, how can the time-variance be formally proven?


texts: Spectral Analysis vs. Spectral Line Analysis

texts:
We usually talk about "spectral analysis" but some resources (this paper or this doc) talk about "spectral line analysis".

Does this make sense to you, i.e. are the 2 fields actually different or the 2 names refer to the same thing?


texts: Calculate input impedance from S-parameters

texts:
I wrote this simple matlab code to calculate input impedance of a low pass filter but not sure whether all the steps are correct. * First attempt

z0           = 50;

A            = importdata("sparameters_lowpass.dat");

freq         = A(:,1);
realPart     = A(:,2);
imagPart     = A(:,3);

s11_freq     = realPart + imagPart*i; % In frequency domain

%X           = fft(s11_freq, N); % Not sure whether this is the right way to interpolate

X            = s11_freq;

X_conj       = conj(flip(X));

X_input      = [X_conj(1:end-1)', X(1), X(2:end)']';

M            = length(X_input);

window       = hamming(M);

s11_time     = ifft(X_input.*window, 1024);

z_in         = z0*(1+ (s11_time))./(1- (s11_time));


texts: Is the expectation of a random process $X(t)$ with zero DC component necessarily zero? Or can it be non-zero depending upon the process?


texts: Have you seen these kind of DFT artifacts?

texts:
we've been working with the Spiral DFT Implementation on an FPGA and managed to get it working.
Unfortunately we see these strange artifacts and have no idea where they come from. Our input is a perfect sine wave, generated by a python script and fed into the simulation as a text file.

We have 2 issues:

Additional peaks in the frequency spectrum that should not be there.
The main amplitude breaks down in the center of the peak.

--> we thought this is due to a value overflow but the phenomenon persists even if we decrease the amplitude of the input signal.

We should note that our DFT takes 16 bit values, which are the upper 16 bit of our 24 bit sample values.

Here are some examples:

Different Frequencies:

Frequency: 10kHz, Amplitude: (2**23-1)/14

20kHz, Amplitude: (2**23-1)/14

Frequency: 30kHz, Amplitude: (2**23-1)/14

Changing the Amplitude results in following:

Frequency: 6.5kHz, Amplitude: (2**23-1)/10

Frequency: 6.5kHz, Amplitude: (2**23-1)/13

Our Sample Rate is 65.1 kHz.

If you have seen this phenomenon before or have suggestions on how to approach this issue, we're happy about any ideas! :)


texts: What happens to frequency domain components if I circular-shift the time domain samples?

texts:
If I have a time domain array of N samples, and swap the final half with initial half, then perform an FFT, what will be the effects on the frequency domain? Will the frequencies' magnitudes remain the same as the ones in the frequency domain of the original array?

What if instead of swapping halves, I just circle 1/10th, or 1/5th, of the full array?

I want to know this because I am trying to take the frequencies' magnitudes at real time, and as new sample data arrives, I thought of simply updating my frequency domain data (subtracting the old samples and adding the new ones) instead of redoing the whole FFT every time.


texts:
I have an image of a cross-section of soil where the main object of interest are the plant roots.

I would like to segment or extract the roots from the background noise of the soil itself and other obstructions such as spider webs. The goal is to form a clean image of the root network or end up with a result similar to what is shown here. In some cases, cracks in the soil itself are detected as interesting lines (structural artifacts) which is an unwanted result.
It seems I may have to apply a series of filters, does anyone have a suggestion what I could do to clean such image(s) or a pre-processing approach that would provide better results.
Current results


texts:
Can anyone please explain about Cepstral Mean Normalization, how the equivalence property of convolution affect this? Is it must to do CMN in MFCC Based Speaker Recognition? Why the property of convolution is the fundamental need for MFCC?

I am very new to this signal processing. Please help


texts: Goertzel filter for non-equidistant sampling

texts:
We have a sensor that reports a new value asynchronously only when the value changes. Therefore, we receive samples at arbitrary time stamps. We basically receive value + timestamp of the measurement at irregular times.
Currently, we re-sample this signal to get equidistant sampling and use the goertzel algorithm to extract a single frequency bin. However, we were wondering if the the Goertzel filter could be extended to work directly on this kind of data.
Does anyone have any experience with this?


texts: Super basic questions on statistical process

texts:
Before starting: I am really a beginner in statistical process in time. I mainly do quantum information and while learning aspect of quantum noise I realized that I am actually too weak on basics of statistical process.

Thus, I would like to first tell you what I basically understand and that you confirm or correct me. Then I have some more specific questions.

Allright, here is what (I think) I understood, and I would like to check

Let's consider a random process. I consider a random variable $X$.
First, I can fix a given time $t$, and ask questions about the statistical properties of $X(t)$. For example its mean value $<X(t)>$, its standard deviation $<(X(t)-<X(t)>)^2>$, and I could imagine many other quantities. In principle without further assumptions I need a lot of information (an infinite one ?) to specify all the characteristics of the statistics of $X(t)$.

Furthermore: here it is what is happening at a given time $t$. Even if I fully described the statistics of $X(t)$ I would lack information. For example there might have some correlations between what happens at $t$ and what happens at $t'$. Thus I also need to specify some properties linking those different times. They are not given by what happens at a given time $t$.

An example I have in mind is a variable that verifies: $X(t) \in [0,1]$, uniformously distributed. But the correlations are such that $X(t'>t)=X(t)$. Then if I do a first experiment I will have $X(t)=0.2$ for example. And for all further times $X(t')=0.2$ as well because of correlations. If I re-run the experiment I can have $X(t)=0.87$ for example but again the values for further times will be the same as $X(t)$. Do you confirm my understanding ?

To specify what happens in time we thus need other information like autocorrelation functions: $C(t,t')=<X(t)X(t')>$. And probably many other quantities that I am not aware of.

More specific questions

What is a rigorous definition of stationnary statistical process ? What I basically understand is that the statistical properties of $X(t)$ actually do not depend on $t$. And the statistical properties of $X(t)X(t')$ only depend on $t-t'$. But I guess there is a more precise definition than this.

In a paper I am reading, they assume a gaussian, stationnary process. They say that knowing the average value $<X>$ and the autocorrelation function allows to fully specify the process. I do not understand this. I agree that it will fully describe the statistical behavior at a given time (because from autocorrelation function we have access to the variance in addition to the mean we already knew). But how is the autocorrelation function the only information we need to specify the statistical relationship between different times ?

Finally, basic question but... what is the name of the precise field that deals with all this ? It is not statistics (this is way too broad), I tried some keyword with "statistical noise" but I didn't find something good. Related: is there a good reference to really give the basics ? I don't think I need much.


texts:
I have found a problem in applying Laplace Transform to $-e^{-at}u(-t)$
I am doing these steps:

$$ = - \int_{-\infty}^{+\infty} e^{-at}u(-t)  e^{-st}dt$$
$$ = - \int_{-\infty}^{0} e^{-at}  e^{-st}dt$$
$$ = - \int_{-\infty}^{0} e^{-(a+s)t}dt$$
$$ = - [-\frac{1}{a+s} e^{-(a+s)t}]|_{-\infty}^{0}$$
$$ = - [-\frac{1}{a+s} (e^{-(a+s)0}-e^{-(a+s)-\infty})]$$

$$ = - [-\frac{1}{a+s} (1- \infty)]$$

$$ = \infty$$
Can anyone help me why it is showing like that.I check it on internet and all the books are showing the answer is $\frac{1}{s+a}$


texts: Finding events in a signal given its number

texts:
I have a signal with14 events

of 2 types:

type 1: the 4th and the 10th peak
type 2: the rest of the peaks

I would like a function in which the input is the number of events of each type  and the output should be the timing and the type of each event.

Does anyone knows any function that performs this (maybe some machine learning algorith). I use R, but I can also use python.

I have signals with events of different shapes (the plot is just one example). I would like to have a function that works independently of the shape of the events so I do not need to provide information about the specific shapes that I have for a given signal.


texts: why does quantisation noise sound like white noise?

texts:
To record sound and transform it to digital form, audio signals are recorded and its pressure levels are sampled and quantised.

Why does quantisation noise sound like white noise?


texts: I read in several sources that the dominant frequencies of the QRS complex in the ECG signal are between more or less 10-40 Hz. What I don't understand is that normally heart beats occur every other second (for a normal heart rate of 60 Beats per minute) - that is a frequency of 1 Hz. In each heart beat there is one QRS complex. So why is the dominant frequency not 1 Hz but 10-40 Hz ?


texts: Continuous Time Signal and Discrete Time Signal - Connection Between Periodicity and Discretness

texts:
I know that all periodic continuous time signal have discrete spectral representations, but are all discrete spectral representations periodic in continuous time?

Also, can all periodic signals be represented by a fourier series?


texts: Filtering data so that only rising edge is left

texts:
I have data that looks like this:

Sometimes the data has a higher point in the middle of the shallow slope

I want to find a way to filter the data such that it smooths it and leaves the first rising edge as the peak like so

Is a weighted average a good way to achieve this? or would that not do the trick? Are there better methods

Right now i have been using a weight average to smooth it but it leaves me with a peak in the center of the slope instead of towards the rising edge.


texts: Filtering and Differentiating phase-modulated signals

texts:
For a project work, I need to demodulate data from a Laser Doppler Vibrometer. The distance information is phase-modulated, therefore the velocity information is frequency-modulated: $$x \propto \varphi \; \Leftrightarrow \; v = \dot{x} \propto \dot{\varphi} = \omega$$

After demodulation the (velocity) signal still is quite noisy, therefore low pass filtering is useful. I've tried filtering both in the distance domain (before differentiating) and in the velocity domain (after differentiating the distance signal). The results seem to be quite similar. What's the difference, exactly?

Specifically, I'd like to know:

What's the unit of the velocity/distance spectrum, exactly?
What is the relation between the (maximum) velocity and the filter's cutoff frequency? Could I specify a "cutoff velocity" instead?
What method is to prefer? Filtering in distance or velocity domain?


texts: Does direct-digital RF signal generation have a "baseband?"

texts: Semiconductor technology has advanced to the point that the digital representation of a modulated RF signal can be applied directly to a DAC for transmission. Is there, then, any distinction between a "baseband" signal and the RF transmitted signal?


texts: Richardson-Lucy iterative deconvolution is producing erroneous deconvoluted images

texts:
I am working on producing deconvoluted image from degraded image where the process of degradation is known to be linear and space invariant. I have used Richardson-Lucy iterative deconvolution method for that. The formula used is given here under.$${\hat{f}}^{k+1}=\left ( {h^{*}}\otimes \frac{g}{h\otimes {\hat{f}}^{k}} \right ). {\hat{f}}^{k}$$ where g is the degraded image, h is the point spread function (PSF) and f's are deconvoluted images. I have tried two processes of degradations (both linear and space invariant). In both cases the RL deconvolution is not working, producing erroneuos images. Although same RL deconvolution script is used in other cases with success. On top of that when I tried Weiner filtering with apodization then in both cases I was able to produce physically meaningful image.
Degradation Type-1
Original activity

Point spread function (PSF)

Degraded image

Richardson-Lucy deconvoluted image iterates (1,3,6,9,12,15,18,21,24,27th iterate images)

Degradation Type-2
Point spread function (PSF)

Degraded image

Richardson-Lucy deconvoluted image iterates (1,3,6,9,12,15,18,21,24,27th iterate images)

Weiner filtering with apodization (for Degradation Type-1)
Original activity image + degraded image + deconvoluted image

The fact that Weiner filtering works put us in a big question of why RL algorithm is not working specially in the scenario where same RL algorithm script is used for successful image production in a case where degradation is less (not shown here). Point to be noticed that in RL algorithm staring from first iteration image it gets worse than degraded image. All the intensity gets accumulated in center region. Why is it happening?
I will be higly obliged if someone provides me hint about where things are wrong. Thanks in Advance.


texts: Common Use Cases for 2D Non Separable Convolution Filters

texts:
In the image processing world, I've noticed that a lot of the popular convolution filters are separable. Here's a quick list of common separable filters:

Sobel
Gaussian blur
Box filter (all ones, for example)

Not separable:

Laplacian
What other common filters are not separable?

I'm a little bit worried about this post being closed for being too broad or not constructive enough. But, some "what is this good for?" questions have been well-received here.

Background: I've developed some high-performance 2D convolution implementations. I'd like to have a better handle on the motivating applications.


texts: RF mask (power spectrum measurement bandwidth)

texts: I'm reading an IEEE standard and it says that the first channel needs to be 50dB below the highest power in a 6 MHz operating channel in 100 KHz bandwidth. It was mentioned several times that the power spectrum measurement needs to be done over a 100 KHz bandwidth. Isn't it sufficient to just say that the power of the first channel should be 50dB below the highest power. I guess my question is what is the significance of the measurement bandwidth, such as the 100 KHz. Why do we need to specify a measurement bandwidth?.


texts: Matching 2 undirected weighted graph in matlab

texts:
Consider 2 undirected weighted graph as shown in figure.Consider first graph as G1 and second graph as G2.Graph G1 consist of vertices V={F1,F2,F3,F4,F5} and Graph G2 consist of vertices V={N1,N2,N3,N4,N5}.

According to logic of my project we have to match vertices of 2 graphs as shown below.

F1=N2
F2=N5
F3=N4
F4=N3
F5=N1

These matching values are calculated manually by seeing values of diagonal elements of 2 adjacency matrices given below.

Adjacency matrix for 2 graphs are shown below.

     F1     F2    F3    F4    F5
 F1   3     4     3     1     4
 F2   3     5     3     1     4
 F3   4     3     3     1     4
 F4   3     2     3     1     4
 F5   3     4     3     1     5

      N1    N2    N3    N4    N5
 N1   4     3     2     1     4
 N2   4     3     3     1     2
 N3   4     3     2     1     2
 N4   4     3     2     1     1
 N5   4     3     2     1     5

How it can be done in matlab?


texts:
Why are those two systems equivalent?

System 1:$$\boxed{x[n]}\longrightarrow\boxed{\uparrow 3}\longrightarrow\boxed{\downarrow 4} \longrightarrow\boxed{y[n]}$$
System 2:$$\boxed{x[n]}\longrightarrow\boxed{\downarrow 4}\longrightarrow\boxed{\uparrow 3} \longrightarrow\boxed{y[n]}$$

By drawing in time the first system's output my result for the first $4$ samples: $$y[0]=x[0], y[1]=0,y[2]=0,y[3]=x[12]$$

by drawing in time the second system's output my result for the first $4$ samples: $$y[0]=x[0], y[1]=0,y[2]=0,y[3]=x[4]$$

Why the statement is correct?


texts: Cross-correlation and auto-correlation of an image using $\tt fft$

texts:
I'm working on a project where I have to find the auto-correlation and cross- correlation of types of pixels in an image:

I am using MATLAB and my code for the autocorrelation is:

fftimage = fft2(image);
congfft = conj(fftimage);%conjugate
ans1 = fftimage*congfft;
ans2 = ifft2(ans1);
ans3 = abs(ans2);
ans3 = fftshift(ans3);
image(ans3,'CDataMapping','scaled')

for the cross-correlation

 ans1 = fftimage*congfft2; %congfft2: conjugate the other type of pixel.

I don't understand clearly the outputs of the code. Here is the autocorrelation output of the yellow pixels:


texts: Map the Amplitude Envelope from [dbmin,dbmax] to [-inf,dbmax]

texts:
I would like to map the envelope of a sound (which has max dB level = dbmax, and min dB level = dbmin)   from [dbmin,dbmax] to [-inf,dbmax].
Here is an example :

It's like some sort of dynamic compression (or rather dynamic range 'expansion' to be precise), but here the mapping is linear from [dbmin,dbmax] to [-inf,dbmax].

What's the simplest way to do this with an  x[n] array? (my input is a soundfile, not an envelope. Is there a simple method for working on envelope (envelope follower ?))

Thank you.

PS : not realtime processing. Thus I can pre-compute dbmax, dbmin.


texts: Is Oversampling a Signal Same as Discretizing the Signal?

texts:
I am confused about the new terminology 'oversampling' I have encountered very recently.
To explain it better, I am going to use numerical examples.
Now let's assume I am using a BPSK modulator and my symbol period is 1 second. To digitize and process this using MATLAB or Python, what I was doing was to discretize the signal with 100. (Using np.repeat() or repmat to information bit sequence).
My questions are:

1) Is this value (100) the oversampling factor?

2) Is this discretization value of mine, is the ratio between symbol duration times sampling frequency of the system ?

3) If the ratio that I mentioned earlier is not what oversampling is, can you explain or give me links for some kind of explanation? ( I found a lot about it but it made me even more confused due to the fact people are using upsampling and oversampling interchangeably, which I think is totally wrong)

Thanks in advance.


texts: Identify fundamental frequencies in short-period waveform

texts:
Maybe I'm being stupid and failing to see something here...

I need to analyse a signal over a relatively short time period (~1-2 cycles for the lowest frequency, perhaps 20-30 cycles for the highest frequency). My aim is to identify the frequencies which are present - usually two or three, never more than six. So naturally, I perform a Fourier transform, applying a Blackman-Harris window to get a reasonable compromise between peak breadth and sidelobe suppression.

A priori I expected that when there were two frequencies I'd see two peaks, with three frequencies I'd see three peaks etc. What actually happens is that with two frequencies I get three peaks ($f_{1}$, $f_{2}$, $f_{1}+f_{2}$), with three frequencies I get up to seven peaks ($f_{1}$, $f_{2}$, $f_{3}$, $f_{1}+f_{2}$, $f_{2}+f_{3}$, $f_{1}+f_{3}$, $f_{1}+f_{2}+f_{3}$) and so on up to 63 peaks with six frequencies. This isn't terribly useful. Because the peaks are relatively broad because of the short sampling time, it is almost always the case that some of them overlap and are indistinguishable from one another, and even if (in the six-frequency case) I could see all 63 peaks, there are $ ^{63}C_{6} = 67,945,521$ possible ways of guessing which six peaks represent the actual frequencies I want!

Yet I can't imagine that what I'm trying to do is that unusual, so I can't help thinking that I'm either doing the transform wrong, or missing some other, much simpler way of achieving my objective. Which is it?


texts: How to filter very low amplitude staircase signal from a background noise

texts:
I have a staircase signal of a specific pattern that is drowned in a background noise. The ratio of signal to noise amplitude is nearly $10^{-4}$. I have attached the image of both with this post. I am an absolute noob in this field so please suggest me a direction on which path I should start to hunt for the solution.

I have tried to do fft on both of the signal and looked for any distinct signature to filter the background noise. Also I have tried with matched filtering and basic convolution techniques. But none of these seem to work. Your advice will be much appreciated. Thanks.


texts: I often see here formula expressed in term of $z$. But what is $z$?

texts:
While searching resources for generating pink noise (and with your help in the comments and answers of other questions), I came to such kind of formula:

  $$ H(z) = { .041 - .096z^{-1} + .051z^{-2} - .004z^{-3} \over{}  -2.495z^{-1} + 2.017z^{-2} - .522z^{-3}} $$

  (from https://dsp.stackexchange.com/a/376/46389)

What is $z$?

I know it should be obvious for anyone familiar with digital signal processing but for people coming from a different field, it is not obvious to know what this represents and how it relates to the signal amplitude or frequency or any other "tangible" parameter.

Edit: amusingly enough, the site suggested me the z-tranform tag for this question. I suspect I should consider that as a clue.


texts: Is there a algorithm for detecting bursts in signal?

texts:
Suppose we have a signal like that and we want to know when a burst have happened. Rigth now our approach is pretty straightforward. We measure std of signal at some time window and then with each new signal bin we measure it's difference with respect to previous and if it's much more than std, then we have found our burst. But maybe this is known problem and there is a library for that, or the way to do it more efficient.

signal time plot


texts:
I am trying to test the sub-pixel accuracy of an optical flow method and need to generate some synthetic data to do so. I start with an image $I_0$, which is a Perlin noise image and a set of sub-pixel ground truth flow fields to apply, e.g. rotation around the image centre. I want to produce the most faithful image, ${I_1}$, that agrees with this sub-pixel flow.

Method 1.
The dumb approach is to just apply the flow directly to the image, and interpolate a new image from the old using one of the many interpolation techniques. But the interpolation is almost certain to introduce noise, the nature of which I guess depends on the statistics of the image and the interpolation function chosen.

Method 2.
A better approach (I think), is to start with a much larger image $J_0$, which is the parent of $I_0$ in that downsampling $J_0$ results in $I_0$. Then apply the flow to the large image and downsample that to get to $I_1$.

However, the downsampling process itself requires an assumption about which kernel to use. Also, this gets tricky if $I_1$ needs to be large (say $1000 \times 1000$), since $J_0$ will need to be even larger (e.g. $20000 \times 20000$) and becomes difficult to store in memory.

Method 3. Another possibility is to model the image as some known continuous function, apply a flow function analytically to the image function, then produce the warped image by sampling from that new function with an assumed point spread function to produce the warped image. But I've no idea if that's possible or what kind of functions would work.

Would anyone be able to provide any suggestions? Would Method 3 work? Is there another, better alternative?

Thanks a lot in advance!


texts: Image processing - Is mean filter shift-invariant?

texts:
Update:
So, according to this article.
For a 2D mean filter,
$$
H[f(x,y)]=\frac{1}{NM}\sum_{k=0}^{N-1}\sum_{p=0}^{M-1}f(x-k,y-p)=g(x,y)\tag{1}
$$
$$
g(x+x_0,y+y_0)=\frac{1}{NM}\sum_{k=0}^{N-1}\sum_{p=0}^{M-1}f(x+x_0-k,y+y_0-p)\tag{2}
$$
$$
H[f(x+x_0,y+y_0)]=\frac{1}{NM}\sum_{k=0}^{N-1}\sum_{p=0}^{M-1}f(x+x_0-k,y+y_0-p)\tag{3}
$$
because (2) equals to (3), a 2D mean filter is shift-invariant.
Is the proof right?
BTW, if f(x,y) is a 2D image, such as

1 2 1
2 3 2
1 2 1

and g(x,y) is (using zero padding)

 8/9 11/9  8/9
11/9 15/9 11/9
 8/9 11/9  8/9

What does these
$$
H[f(x+x_0,y+y_0)]\tag{4}
$$
and
$$
g(x+x_0,y+y_0)\tag{5}
$$
mean?
(I want to prove it by using a image.)
Origin:
I read this article
and this Wikipedia page.
However, I still don't know what it means for (or how to apply to) a mean filter.
Can anyone give me an image example?
(for example, https://dsp.stackexchange.com/a/14435 is easier to understand)


texts: Windowing in frequency domain. how to cut a spectrum using window function in frequency domain?

texts:
I have a set of measurement data of a FMCW radar and trying to estimate the system transfer function. In the frequency spectrum, I can see multiple peaks due to the non-linearity that exist in the system. I want to estimate the transfer function for each peak.

What I did:
I designed filters to cut out each peak and using deconvolution I've estimated the transfer function. In the verification step somehow this process has failed to replicate the system exactly.

What I wanted to do:
I want to directly cut out frequency peaks using window function instead of filters. But I am not able to understand how to place the window around frequency peaks so that I can extract each peak separately.


texts: How to calculate the impulse response from the transmission function of an optical device?

texts:
Let $f(x,y)$ be the electric field in complex representation processed by an optical device described by a transmission function $t(x,y)$ such as:
$$
f'(x,y)=f(x,y)t(x,y)
$$

Where $f'$ is the output electric field.

How can I calculate the impulse response of the device? I don't understand how to determine a function $h(x,y)$ such as:
$$
f'(x,y) = f(x,y)*h(x,y)
$$

More specifically I'm working with a lens transmission function in the paraxial approximation:
$$
t(x,y)= P(x,y)\exp\left[-\frac{jk}{2f}(x^2+y^2)\right]
$$
where $P(x,y)$ is the pupil function of the lens (usually a $circ(x,y)$ function).

P.S.
I wonder if there is a general solution to the mathematical problem of finding a function $h$ for a given $t$ which satisfy $ft=f*h$ for a generic function $f$.

I've posted the same question here: https://physics.stackexchange.com/questions/304444/ but I haven't received any answer.


texts:
I want to prove using numpy the theory of Fourier transforms in which translation in space corresponds to a shift in the phase domain (frequency domain remains constant).
I have generated these three 2d-numpy arrays showing translated versions of the same boxes:

I then use numpy to calculate the FFT of each of these images doing:
H1 = np.fft.fft2(im1)
H2 = np.fft.fft2(im2)
H3 = np.fft.fft2(im3)

Although the frequency and phase plots of the transformed arrays look very similar I perform a comparison (MSE) between them to see the differences.
This is the MSE plot for the frequency domain:

And this is the MSE plot for the phase domain:

I'm confused as I was expecting to see identical spectra in the frequency domain but different in the phase domain.
I would appreciate if someone can provide some insight to detect if the problem comes from a theoretical misunderstanding or a problem in the code.
Here is the code to replicate the results shown here.


texts: Higher-order splines all mostly look like Gaussians. I wonder how higher-order splines (from 3rd order onwards) and Gaussians differ.  How can one  differentiate them?


texts:
I am still new to MATLAB, so apologies if I sound lazy to you. I am attempting to model a transformation as a set of matrix operations. I start with a vector, up-sample it by $U$ (up-sampling rate), take its Fourier transform, apply a Low-pass filter, take Inverse Fourier transform and finally down-sample at a rate $D$, where $U \neq D$.

I want to do all these transformations as a set of matrix multiplication, and not use MATLAB's built-in commands. Up-sampling / Down-sampling is easy. And so is Low Pass Filtering. But I am stuck with fft and ifft. I could use dftmtx to generate the matrix that would perform the fft, but using a vector $1000$ samples long, and up-sampling it by $1000$ makes a very big array, so that matrix generated by dftmtx becomes memory constrained. I have looked the code for dftmtx which is simply fft(eye(n)), where $n$ is length of the vector to be transformed. I tried replacing eye(n) with speye(n), so as to reduce the matrix values, but it would not work this way too.

The second problem is I do not have any idea how to do ifft using a matrix multiplication. I understand this matrix will also be a very large matrix and become memory constrained. However, no search on any online resource has helped me so far.

Whereas doing the above mentioned transformations is easily done by MATLAB, I need the matrix version for my application.

Any guidance is appreciated.

My code in simplest for is

n=1000;% length of original code
a=randn(n,1)+j*randn(n,1); %initialise
U=1001; %up-sampling rate
u=zeros(n*U,n); %Blank up-sample matrix
%find indices for up-sampling
for i=1:N;
    u(U*i-(U-1),i)=1;
end
ua=u*a; %up-sample a by rate of U
%Find Fourier transform
ft=dftmtx(length(ua));
fa=ft*ua;
%Designing LPF
lp=speye(length(fa)); %initialise
j=length(fa)/2-ceil(U/3):length(fa)/2+ceil(U/3); %define central region where high frequencies will be ignored
lp(j,j)=0; %ignore center values which have high frequency, render this region zero
lpa= lp*fa; %Low Pass Filter Up-sampled Fourier transformed sequence
%Inverse Fourier transform
ifa=ifft(lpa); %How to implement this as matrix multiplication so that we have ifa = ifft_matrix * lpa
%Down-sample
D=1000; %down-sample rate
j=length(ifa)/D;
d=zeros(j,length(ifa)); %Blank up-sample matrix %find indices for down-sampling code
for i=1:j
    d(i,D*i-(D-1))=1;
end
da=d*ifa; %calculate the down-sampled matrix da with length 1001 (original code stretched by 1 sample)


texts: Python butter filter: choosing between analog and digital filter types

texts:
I have a system that performs wireless sampling of data about every 7.5ms (133Hz). Due to it being wireless, I get occasional data drop out. I want to construct a LP butter filter with cut-off frequency of 10Hz using Python's scipy butter method and then downsample everything to a lower frequency.

One of the options that it asks for is whether I want to use analog or digital filter. Isn't running this filter offline in python automatically assume that it's digital? Or does the fact that my samples don't necessarily come in at regular intervals mean that I should stick with analog filter?

I understand how a physical analog filter differs from a digital one, but does this difference apply in the same way to a python-based butter filter?


texts: fundamental question regarding bandwidth

texts:
What I've learnt so far is that, if we change a time-domain signal into frequency domain, its bandwidth shows the range of frequency of radio waves that required to generate such signal.
However, for example in fiber optic,a 1550nm light is used to be the source and BPSK is used to modulate it, but if there is only a single frequency of light (1550nm) forming a signal, then what is the point of finding the bandwidth?


texts: Effect on T1 relaxation due to the application of the gradient

texts:
In MRI, the T1 relaxation constant for a tissue type is defined as the time taken for the longitudinal relaxation to recover approximately to 63% of its initial value.

Is there any modification of the T1 time due to application of gradients. Suppose, after application of the RF pulse, the magnetization is tipped completely in the transverse place. If one waits, then 63% of the longitudinal magnetization (z-direction) will recover after T1 seconds. Now, say at the same time I apply a gradient along the x or y direction. Would this have any affect on the longitudinal relaxation time?


texts:
Let's assume I have the following transfer function:

$$
H(z)=\frac{z-\left(\frac{1}{\sqrt{2}}+i \cdot \frac{1}{\sqrt{2}}\right)}{z}
$$

It looks like a first order highpass-filter with a complex zero at $\frac{1}{\sqrt{2}} + i \cdot \frac{1}{\sqrt{2}}$.

If I do the inverse $z$-transform, I get the following difference equation:

$$
y(k)=x(k) - \frac{1}{\sqrt{2}} \cdot x(k-1) - i \cdot \frac{1}{\sqrt{2}} \cdot x(k-1)
$$
with $k$ as time step variable. My input signal is real valued and my output should be real valued as well.

How do I implement that filter e.g. in C?
How do I deal with that complex numbers in that case?


texts:
I perform FFT of two real sequences using "two for the price of one" - array $x[n]$ and $y[n]$ is packed as a complex number - $z[n] = x[n] + j y[n]$ - using FFT2 and FFTW. I am able to decode the half-complex output format of FFTW and the result match with the output of FFT2.

The output conversion goes like this. If FFT of complex sequence $z[n]$ is $Z[k]$ where $Z[k] = FFT\{z[n]\}$. The equations look something like this

$$FFT\{x[n]\} = (Z[k] + Z^*[N-k]) / 2 $$

$$FFT\{y[n]\} = (Z[k] - Z^*[N-k]) / (2j) $$

When I perform inverse FFT with FFTW I am unable to understand the output format of FFTW and the inverse output or FFT2 do not agree with FFTW. Any clue on   how to transform the output of inverse FFT to generate inverse FFT of two sequences.


texts: Channel Vocoder producing output with "click" sounds

texts:
I'm trying to make a channel vocoder that takes two inputs, one a frequency rich carrier (a musical sound) and the other a modulator (vocals). Applying operations involved in the channel vocoder generates an output that makes the vocals sing in the tone of the carrier. (Think Songify.)

Presently I'm splitting both the signals into windows of a specified size and also overlapping two consecutive windows by a certain factor. FFTs of the windows of the carrier and modulator are then computed. The fft is further split into bands of equal length. The average of magnitudes of modulator's fft in that particular

for band_no = 0 to band_count:
        for i = 0 to band_length:
                avg+=magnitudeOfModulatorFFT[band_no*band_length + i];

I am then multiplying "avg" to each of carrierFFT's frequency value lying within the band, and then taking the inverse FFT. Although this produces the desired output, it also includes a clicking sound. Even after multiplying both the signals by a window function.

Any way I could solve this?


texts: how we can prove that the Variance of (white Gussian noise) is average power of the noise?

texts: White Gaussian Noise PSD is fixed No/2. So how we can investigate that noise variance is average power of the noise?


texts: Generic response of an IIR filter from its transfer function

texts:
How can I derive the general response of an IIR filter from its transfer function? I know that:

$$H(z)=\frac{1}{1+\sum\limits_{m=1}^N{a_m z^{-m}}}$$

Thus:

$$Y(z)=X(z)H(z)=\frac{X(z)}{1+\sum\limits_{m=1}^N{a_m z^{-m}}}$$

The general response is:

$$y[n]=-\sum\limits_{m=1}^N{a_m y[n-m]}+x[n]$$

Where the above expression comes from?

Thank you for your time.


texts:
The Wikipedia page on sound localisation states:

  Localization accuracy is 1 degree for sources in front of the listener
  and 15 degrees for sources to the sides. Humans can discern interaural
  time differences of 10 microseconds or less.

However, I can't find any more detailed information on localisation accuracy.

What I'm after is a plot of accuracy (Y) against azimuth (X).

Can anyone get me one step towards this information?

π


texts:
Since we know about convolution filter in image processing, it is basically a matrix applied to an image and a mathematical operation comprised of integers and it works by determining the value of a central pixel by adding the weighted values of all its neighbours together.
Question: In Envi software, we apply a convolution filter by going under the Filter  >  Convolutions and Morphologytab, Now there we have a default 3×3 LPF (low-pass filter) and HPF(High pass filter), On the basis of their usage/kernels, how exactly these filters works?
Note: I am a beginner in image processing so need a basic explanation.


texts:
Given a set of $N$ iid samples, $x_k$, of random variable $X$ and the expectation estimator,
$$E[X] \approx \frac{1}{N} \sum_{k=1}^N x_k,$$
how do I find the variance of this estimator (not the variance of $X$) as a function of $N$?
I would only like to know the error in the estimate of $E[X]$ and what I would expect is that more samples (larger $N$) drives the variance down, so that
$$\lim_{N \to \infty} \bigg[\frac{1}{N} \sum_{k=1}^N x_k - E[X]\bigg] = 0.$$
Am I right about this? Perhaps variance is not the correct way to measure the error in the estimate? What I would like to understand is, how many $N$ is sufficient because there is a cost for every extra $x_k$ that is processed.


texts: What is the physical interpretation of the dB scale on a bode plot and what is a negative dB?

texts: I have no physical interpretation of the Bode plot. What does it mean for a bode plot to have negative dB over its entire duration on the log-scale frequency?


texts: Libraries for Voice Activity Detection (Not Speech Recognition)

texts:
As follow up to my previous question I was wondering if there are any speech detection libraries in existence. By speech detection I mean passing in an audio buffer and getting back an index of where speech starts and stops. So if I have 10 seconds of audio sampling at 44kHz, I would expect an array of numbers such as:

44000
88000
123000
190334
...

This would indicate for example that speech starts one second in and then finishes at the two second point, etc.

What I'm not looking for is speech recognition which writes out text from spoken word. This unfortunately is what I see a lot of when I google 'speech detection'.

It would be great if the library was in C, C++ or even Objective-C as I'm writing an app for the iPhone.

Thanks!


texts: I tried to find the stopband edge frequency.I know that the signal amplitude falls to -3 dB below the passband edge frequency.I examine a lot of DSP books so i don't any information about how to calculate stopband edge frequency.


texts: Low-pass filter for cleaning vocal speech before performing onset detection

texts:
I am creating a piece of software that detects when speech starts in an audio recording. Before performing onset detection I would like to remove any noise that may cause any false positives. Which would be the best low-pass filter to use for this?

EDIT 1

Maybe I went too far in saying that I am performing onset detection. Basically all I am doing is looking for the first time the amplitude goes above a certain value and taking that as a flag for the start of the envelope. It is simple and crude, but I believe if my low-pass filter successfully removes noise then the onset detection should work. Or am I totally off on this?


texts: Show Equivalence Between Multiplication in Time Domain to Convolution in Frequency Domain

texts:
My goal is to compute the Fourier of the product between two discrete time signals, y1 and y2. This can be done by computing the convolution between the fourier transform of y1, f1 and the fourier transform of y2, f2 (This is what I understood from the Wikipedia Page).

I tried the following code in MATLAB, it shows that the convolution between y1 and y2 is equal to ifft(f1.*f2) as shown on the Wikipedia page, the code worked perfectly:

y1=sin(x);
y2=sin(3*x);
n=length(y1);
convol=conv(y1,y2);
f1=fft(y1,n*2-1);
f2=fft(y2,n*2-1);
plot(ifft(f1.*f2))
hold on
plot(convol,'x')

Now I want to be able to compute the fourier transform of the product of y1 by y2 using the convolution theorem, this can be done by computing the convolution between f1 and f2, so I wrote the following code:

x=0:0.01:10;
y1=sin(x);
y2=sin(3*x);
n=length(y1);
f1=fft(y1,n*2-1);
f2=fft(y2,n*2-1);
convol=conv(f1,f2);
ffty1y2=fft(y1.*y2,length(convol));
plot(abs(ffty1y2))
hold on
plot(abs(convol)/length(convol),'x')

As you can see the results do not match (neither in location of the peaks nor in their values), what could be wrong here? How do I fix this?


texts: How does a FIR equiripple filter behave close to DC?

texts:
I am multiplying two sine waves with the same frequency (f), but might have a phase difference 0 <-> 90 degrees. The product is a two frequency sinusoidal wave with f1 = f-f = DC and f2 = f+f = 2f. I now want to filter out the 2f component. I am currently sampling at Fs = 32 x f. I get a 40th order equiripple FIR filter in Matlab which is ok, but if I later am to increase the sampling rate to 64xf, 128xf etc the filter order gets very high.

I assume however that f1 will not be 0 Hz exactly but somewhere close to 0 Hz. Over this region the ripple in the passband might be unwanted. Can the ripple change that much over such a narrow range? (within whatever limit one has set as the maximum ripple in the passband).

I am using the Filter Design & Analysis tool in Matlab and saw that it does not seem like an equiripple filter has any ripple close to DC (0<->1Hz)? Is this correct?

I need a filter that can attenuate the 2f component with 60dB and need as little computing power as possible. (working on a MCU without FPU using fixed-point).

Is there any other filter that can do a better job than equiripple for what I described? I am open for other types of filters as IIR if someone can convince me of their advantages.


texts: From my knowledge of stability, I understand that if the function approaches a finite number then the system will be stable. Thus if a system is stable its poles will be on the left of the $j\omega$-axis. This means that $Re(p_{1}, p_{2}) < 0 $. What I don't understand is how you determine whether the Imaginary part is equal to 0 or not.


texts:
Usually, algorithms as SIFT, SURF and many others provides a set of $k$ keypoints and the associated descriptor in $d$ dimensions (for example, in SIFT each descriptor has $d$=128 dimensions).

So, in order to describe an image we need a matrix $k \times d$ ($k$ descriptor vectors, each one in $d$ dimensions). So far so good.

My question is: how can we describe an image through a single vector?

This could be really useful since we could save a lot of space and because certain algorithms (like LSH) requires a vector as input/query.

In some papers (for example this, section 6.5) this approach is described as "global descriptors".

Up to know, I found only this paper but it doesn't seem so accurate (and it's from 2009, not so new).

APPLICATION:

I'm trying to implement a memoization optimization to the following application: given an image (taken with a smartphone) containing a movie poster, return the internet movie database's score of that movie.

The optimization (so what I'm trying to implement) is about "remembering" (memoize) the input image and the relative result (IMD score) and for each new input image, search "the most similar one" and return the already computed result (avoid the database query) if a certain threshold is reached.

Two "similar images" in this context are two images pictures containing the same poster but from different prospective and/or places.

The actual approach is solving an approximate nearest neighbor through LSH in order to find the most similar image, but in order to do that I need to represent a query (an image) through a vector, and not through a matrix.


texts: What are some suitable methods for removing low-frequency line noise from measured data prior to system identification/deconvolution?

texts:
I am performing a system identification by striking an experimental model with an instrumented impact hammer and measuring the strain response.

The timebase for the recorded signals is 5 microseconds and the hammer pulse width is typically around 200 microseconds. The data are recorded for 20 milliseconds prior to the hammer hit and for 60 milliseconds after.

An example of the recorded data is provided here. A close-up of the pre-trigger strain signal is provided here.

Prior to the hammer hit, the strain signal is non-zero with both low and high frequency noise present. The low frequency noise is of most concern as this affects the quality of the system impulse response function that I subsequently compute.

Unfortunately when completing the tests I did not record a zero-reference prior to each hammer hit. Thus my questions are:

1) What methods are available for removing the low-frequency line noise from the entire signal given that the strain response is also of low-frequency.

2) Is is possible to use the first 20 milliseconds to construct a longer zero-reference signal that I may then subtract from the recorded data? If so, how?


texts:
I have an image with arbitrary area (the boundary is shown in gray):

I would like to blur just the image within boundaries.

Ordinary blur causes the background proliferating to image and vice versa. This results in dark areas near boundaries of the image:


texts: Are there analog re-configurable filters?

texts:
I know that we can build a chip (or circuit) or use FPGA to have re-configurable digital filters. But how about analog filters, can we have re-configurable analog filters?
thanks


texts:
I read in a paper http://tcts.fpms.ac.be/publications/regpapers/2011/nrmmbgtd_icvs_LNCS.pdf about spatio-temporal filters, which is applied to frames of the videos. I could not understand these terms in the paper if anyone can help.

3D Low Pass Filtering

Temporal Neighborhood

Spatio Temporal Filters


texts: How to transmit voice signal and receive on radio?

texts:
I want to transmit data like radio stations do.

Like radio stations broadcast music and etc. using frequency, so I am also curious to know about it that actually how it works.

I also want to make such a device that will broadcast my data or voice in my apartment and if anybody will on that frequency they should hear my broadcast data.

Is it possible or not ?
If yes then anyone please guide me.
And what happen if i ll broadcast on same frequency on which already data is broadcasting ?

I want to know about radio signals, so any guidance on it.
Thank You.


texts: Obtain transfer function from signal flow graph

texts:
I have a pretty complex signal flow graph with all the blocks and the branches completely describing an analog filter, in the Laplace complex variable s. The signal $X(s)$ is the input and $Y(s)$ is the output.

Is there any (free) tool which can analyze such a graph and provide the overall transfer function $H(s) = Y(s)/X(s)$?

Googling I found this Matlab tool, but the tool it is not free; anyway is there anything similar, in Matlab or in other environments?


texts: How to understand the discrete and continuous convolution in matlab?

texts:
For example,I want to calculate the convolution of two continuous signal in matlab as the image shows.
Indeed, I chose to discrete the signal to compute in matlab.

Obiviously, the 1st and 3rd are right, I think the difference of discrete and
continous convolution might be the reason, but why?


texts:
In this answer, an iterative algorithm is given for recovering a signal from its spectrogram, assuming overlapping windows:

$$
x_{n+1} = \operatorname{istft}(S \cdot \exp(i \cdot \operatorname{angle}(\operatorname{stft}(x_n)))),
$$

where $S$ is the spectrogram and $\operatorname{(i)stft}$ is the (inverse) short-time Fourier transform. I tried this out using NumPy, SciPy, and librosa:

import librosa
import numpy as np
import scipy

def mse(phase_pred, phase_true):
    '''
    Calculate the mean square error between the true phase and the
    predicted (reconstructed) phase.
    '''
    return np.mean(np.angle(phase_pred/phase_true)**2)

# Load an audio file and calculate STFT.
x, sample_rate = librosa.load('audio.wav', sr=44100)
D = librosa.stft(x)
mag, actual_phase = librosa.magphase(D)

# Try to reconstruct the phase using the iterative algorithm above.
phase = np.exp(1.j * np.random.uniform(0., 2*np.pi, size=actual_phase.shape))
x_ = librosa.istft(mag * phase)
print('iter {} mse {}'.format(-1, mse(phase, actual_phase)))
for i in range(100+1):
    _, phase = librosa.magphase(librosa.stft(x_))
    x_ = librosa.istft(mag * phase)
    print('iter {} mse {}'.format(i, mse(phase, actual_phase)))
    if i % 10 == 0:
        scipy.io.wavfile.write('recons{:05d}.wav'.format(i), 44100, x_)

As far as I can tell, it isn't converging:

However, the audio after the 100th iteration definitely sounds closer to the original than the one with randomized phases.

Why isn't this algorithm converging to the correct phase? Am I misunderstanding its purpose?


texts: Matrices of complex numbers multiplication

texts:
I'm trying to implement the multiplication of two matrices something like this picture in c langage.
I want to read the numbers from a text file of x and store it later in an array the code that i tried is below
So the problem that I have is when i use the  static method Complex x[4] = { { 1,0 },{ 0,-1 },{ 1,1 },{ 4,-1 } }; it's work but when i want to read it from input.txt Complex x[4] it does not work what i missed.
typedef struct Complex
{
    double real;
    double image;
}Complex;

int main()
{
   Complex x[4] = { { 1,0 },{ 0,-1 },{ 1,1 },{ 4,-1 } };
   Complex y[4] = { { 0,1 },{ 1,-1 },{ 2,-3 },{ 4,0 } };
   FILE *input;
   input = fopen("input.txt", "r");
   if (sensorInput == NULL ) {
        printf(stderr, "Cannot open file\n");
        return 1;
    }
    else {
        for (i = 0; i<4; ++i)
        {
                fscanf(sensorInput, "%lf + j*%lf\t", &x[i]);
                printf("\n");
        }
    }

}

input.txt
1,0
0,-1
1,1
4,-1


texts: How to quantify delay in filtered output?

texts:
I have a signal containing two sine tones 30 kHz and 60 KHz. I used a second order IIR notch filter to filter out the 60 KHz component from the signal, so now the signal has only 30 KHz component.

When comparing the filtered out signal with a pure 30 KHz tone generated, the earlier samples are very much different but later samples are matching to some good extent.

Plot with earlier samples:

Plot with later samples:

What are the reasons why there is a delay before response being proper?
How can i mathematically quantify this?

Like after how many samples the filtered output will be reasonably matching with that of pure tone?


texts:
I have a device with an accelerometer and I want to make it as precise as possible (using algorithms). Currently no DSP algorithm is running - I am reading values as they are. I've made 3 tests - in each I rotated the device around one of the axes and recorded the values coming out from accelerometer. It would be perfect if each of the tests would produce a perfect circle; this would mean no further calibration is needed. The desired precision is +- 10mg (+- 1%). For example, the mean value when Z axis is maximal is 1024mg. 1010mg would be acceptable.

Below are the results of tests (Matlab copy-paste to display). Is it enough to test the three axes or should I also test the cases where the gravity gets distributed in all of the axes? Further; what approach would you suggest considering calibration to 1% precision (please be specific).

I need a thorough answer on this topic.

Rotation around X axis:

Round Y axis:

Round Z axis:


texts: What is the peak resonance of convolving with a sine FIR filter?

texts:
I'm trying to improve my understanding of FIR filters. As an experiment, I've manually created an FIR filter, whose coefficients follow exactly one period of a sine wave. I'm wondering what is the peak gain frequency of such a filter.

For example, I'm filtering a 44.1 kHz audio signal with filter coefficients corresponding to a 440 Hz sine wave1. Intuitively I would expect that the convolution reinforces this frequency. However, I noticed that the frequencies of slightly lower frequencies (roughly a few semitones, ~365 Hz) actually have a stronger gain. This is also confirmed by running the filter coefficients through scipy.signal.freqz (the red line corresponds to 440 Hz, the gray lines indicate frequencies of neighboring semitones):

Most likely I've encountered a basic rule of digital filter design that I'm not aware of. I'm wondering:

Why is the peak shifted w.r.t. the frequency that is used in the convolution?
Is there a formula that describes the relationship between the two, i.e., is it possible to compute the frequency that the sine convolution should have so that the frequency response peaks at 440 Hz?

1 With a sampling frequency of 44.1 kHz, the period of a 440 Hz sine has a fractional length 100.23 samples, which I'm rounding down to 100 filter coefficients. Could this imperfection play a role in the effect I'm seeing?


texts: Autocorrelation of a stochastic process which is a sum

texts:
Let's say that a stochastic process is defined as
$$X(t) = \sum_{n = -\infty}^{+\infty} X_n * f(t- nT)$$
where $X_n$ is the $n$-th symbol among independent random symbols. I have proved that $\mathbb E[X(t)] = 0$.  But, I want to prove that this stochastic process is Wide Sense Stationary(WSS) which means by definition that its autocorrelation function $R_{xx}(t+k,t) = \mathbb E[X(t+k)*X(t)] = R_x(k)$ i.e. the auto-correlation function of $X(t)$ depends only on $k$, which is the time difference between $t+k$ and $k$. What I have proved is that the auto-correlation function of $X(t)$ is zero but I am not sure that this proves that this stochastic process is WSS. I have attached a picture of the problem and what I have done so far in order to be clearer.


texts: pyo: Retrieve amplitudes of audio signal frequenciy ranges

texts:
I'm playing around with audio signal processing for the first time and just to start somewhere, I decided to use the Python library pyo. My goal is to analyze an audio input signal from the microphone to detect different types of strikes on a drum and to emit an event (or call a callback, whatever) in that case. My first step is to define some frequency ranges, I want to observe and then to have a function called periodically and pass in the peak amplitudes within those defined frequency ranges.

I found that there is something like FFT/DFT and my naive idea is, that I can setup an FFT object with my input signal and create an instance of the Pattern class, that calls a function periodically. In that function I can then call the FFT object in a way, so that it will return me a list with 1024 items. The "whole" frequency spectrum is sliced into 1024 slices and each item in the list contains an amplitude value for that specific frequency range slice at the given time.

I tried the following code and now I don't know how to interpretate the values provided by the FFT.get() method. My question is: Are my expectations wrong or how do I use the FFT class to get my list?

from pyo import *

server = Server(nchnls=1, buffersize=128, duplex=1, winhost="asio").boot().start()

class FrequencyRangeAmplitudeProvider:

    __frequency_ranges = ()
    __listener = lambda x: x

    __fft = None
    __pattern = None

    def __init__(self, input, frequency_ranges, listener):
        self.__frequency_ranges = frequency_ranges
        self.__listener = listener

        self.__fft = FFT(input, size=1024, overlaps=4, wintype=0)
        self.__pattern = Pattern(function=self.__extract_and_emit, time=0.06).play()

    def __extract_and_emit(self):

        #print(self.__fft.get('real', True))
        #print(self.__fft.get('imag', True))
        #print(self.__fft.get('bin', True))

        my_list = []  # <-- How to get the data, I'd need?

        peak_amplitudes = self.__extract_peak_amplitude_in_frequency_ranges(my_list, self.__frequency_ranges)
        self.__listener(peak_amplitudes)

    @staticmethod
    def __extract_peak_amplitude_in_frequency_ranges(amplitudes, frequency_ranges):
        return (('A', -2.34), ('B', -.65))

frequency_ranges_to_observe = (('A', 430, 450), ('B', 870, 890))

def do_something_with_peak_amlitudes_in_frequency_ranges(amplitudes):
    print(amplitudes)

provider = FrequencyRangeAmplitudeProvider(
    Input(),
    frequency_ranges_to_observe,
    do_something_with_peak_amlitudes_in_frequency_ranges)

server.gui(locals())
```


texts: Why is LPCC and LPC used in speech recognition

texts:
I am currently trying to understand why LPCC and LPC is being used in speech recognition and why?

LPC as I understand is model similar to source filter model, modelling speech production,
and LPC thereby tries to estimate the given speech signal, with minimal error.

So we having an estimated signal of the original… but how is that useful in speech recognition task?


texts:
Here's an image of a cracked pattern. I wish to convert it to binary image where the continuous crack lines form the only discernible pixels. The main problems I am having in the lack of uniform lightning, which confuses threshold tools in image processing software, and the "black dots" which should be filtered out.

I am trying to use either Gimp, imagemagick or ImageJ to accomplish this, but have so far not been able to get to decent result.


texts: Where do I start in learning how to demodulate a BPSK31 signal in software?

texts:
I'm a software developer and a Amateur Radio operator. Everything I know about audio (which isn't saying much) is self taught.

I've learned how to generate audio samples and build wav files, send those samples to hardware output, etc. I even figured out how to generate BPSK31.

Where I'm stuck now is demodulation. I've read around enough to realize that this is a much harder side of the problem than modulation, and I understand intuitively why this is.

My question is: Where do I even start to learn how to do this? I'm not even sure what the right terms to search for in Google are. So many of the results are hardware solutions or new theoretical approaches to the problem that I've gotten lost in the noise. What should I be looking for to figure this out?


texts:
I've a noisy signal, which is the sound of motor with a constant speed, so  the sound "should" be  periodic, I know that there is a way to use the autocorrelation function to get the period,I did it, but I can't figure out the period. Any idea how to do that bellow the signal and the result of the autocorrelation :
the signal

the autocorrelation result :


texts: How to specify the input data bit width for fixed point implementation in windowed OFDM?

texts:
I am implementing a fixed point model of Windowed OFDM and to quantize the input data, I used 16 bits(complex) to encode it. But I am not sure how do we specify the minimum number of bits that should be used to encode the input ensuring the output wont be "off" from the expected results.

I saw a similar question in Bit width for FFT but the person has said assuming "M" input bits. Is it an empirical data that you set by trial and error or is there some maths behind it.

Thank you


texts: Can we calculate moments (mean, var, skewness, kurtosis) of a signal in frequency domain?

texts:
The computation of the mean, variance, skewness and kurtosis of a time series in the time domain is straightforward from their definition formulae.

I need these statistic values in different frequency bands. Currently I apply different bandpass filters to the signal and repeat to compute the statistics. I'm curious if it is possible to compute them from the spectrum in the frequency domain?


texts: How can I perform object tracking using noisy images?

texts:
The problem

Take 50+ of Lena images, add gaussian noise to them to the point that you couldn't say it is Lena. I insist, gaussian noise, not gaussian blurring.

Now imagine a cropping window randomly moving around Lena's face (like a video camera focused on this part).

Now you have a handful of noisy pictures containing what used to be Lena's face somewhere in the image.

The question

Now I want to get back the face without knowing it in advance :)

First I have to identify overlapping parts (i.e. noisy faces here) from the set of images.
Then I must superimpose the faces.
And finally, add them to get back her beautiful face.

The problem is that I have no clue of how to do that (cross correlation, perhaps to identify the faces and the shift between them).
I've never had courses in image processing. I just need a reference or an algorithm name.


texts: How to find correlation between two signals regarding peaks?

texts:
In a super-resolution task, I try to predict the higher frequencies of a signal.
X low-resolution (lr) is my input and Y high-resolution (hr) is the target. Where X contains only frequencies till 2 kHz and Y till 4 kHz.

I would like to know which frequencies in Y have a correlation with frequencies in X.
How can I find peaks in both areas (lr and hr) which occur always together?
for example, a peak at 1kHz always means there is also a peak at 3kHZ. In general, I would like to find patterns in the low-resolution input that helps me to find a correlation in high-resolution.

My low-resolution signal is generated from the high-resolution signal. It does not contain aliasing, because it was low-pass-filtered and downsampled correctly.

As an example plot : the lr and hr spectrum. Lr is generated by low-pass-filtering the hr spectrum. I would like to predict the higher frequencies.The lr image does not contain aliasing. So i need to find patterns in the lr signal for example maybe peaks, which always occur together with peaks in the hr signal. That would mean that there is a correlation between lr and hr.


texts: What is the correct order of operations for a 2D Haar wavelet decomposition?

texts:
The source code of iqdb contains a 2D Haar transform implementation. The author claims to have implemented it according to the paper "Fast Multiresolution Image Querying", which is freely available here.

This is the relevant text from the paper:

  A standard two-dimensional Haar wavelet decomposition of an image is very simple to code. It involves a one-dimensional decomposition on each row of the image, followed by a one-dimensional decomposition on each column of the result.

  The following pseudocode performs this one-dimensional decomposition on an array A of h elements, with h a power of two:

proc DecomposeArray (A : array[0..h-1] of color):
    A <- A / sqrt(h)
    while h > 1 do:
        h <- h/2
        for i <- 0 to h-1 do:
            A'[i] <- (A[2i] + A[2i + 1]) / sqrt(2)
            A'[h+i] <- (A[2i] - A[2i + 1]) / sqrt(2)
        end for
        A <- A'
    end while
end proc

  In the pseudocode above, the entries of A are assumed to be 3dimensional color components, each in the range [0,1]. The various arithmetic operations are performed on the separate color components individually.

  An entire r x r image T can thus be decomposed as follows:

proc DecomposeImage(T : array[0..r-1, 0..r-1] of color):
    for row <- 1 to r do:
        DecomposeArray(T[row, 0..r-1])
    end for
    for col <- 1 to r do:
        DecomposeArray(T[0..r-1, col])
    end for
end proc

(end quote)

Implementing it this way does not produce results that match the example images in the majority of the articles I have found on the internet covering this topic, including the Wikipedia article.

Note: The image is divided into 4 large squares, and (only) the top left square is further divided into 4 squares.

However, I have also found counterexamples (i.e. examples that follow the scheme used in the paper above), e.g. here.

The question is whether to loop over all rows and columns, and doing the fully recursive transform in an inner loop for each row or column, - OR - doing one pass of recursion in the outermost loop, and within each single pass only processing the remaining rows and columns.

I have implemented both approaches to demonstrate the difference visually: https://bplu4t2f.github.io/wavelet_toy/

In the approach that wikipedia uses (which I call "pass major" because the pass of recursion is the outermost loop), the emerging pattern shows that each pass divides the image into 4 squares, and only the upper left square is modified in the next pass.

In the approach that iqdb uses (which I call "pass minor"), the emerging pattern shows that only the bottom right of the 4 divisions remains unchanged during subsequent passes.

The pass minor approach feels incorrect to me, because, when looking at it intuitively, it reprocesses parts of the image that have already been transformed during each pass, effectively applying a primitive edge detection scheme recursively on previosly detected edges. It doesn't seem to make much sense to me.

Which of these approaches is correctly being referred to as 2D Haar wavelet decomposition? To both approaches have a name?


texts: How is compression ratio related to number of dct coefficients

texts:
Lets say we have a picture with 6bits/pixel and we use discrete cosine transform with 8x8 blocks of pixels and each coefficient is 8 bits. We can use 2 8x8 masks, a zonal mask A and a zonal bit allocation mask B.

A = [1 1 1 1 1 0 0 0 1 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...... ]
B = [8 7 6 4 3 2 1 0 7 6 5 4 3 2 1 0 6 5 4 3 3 1 1 0 4 4 3 3 2 1 0 0 3 3 3 2 1 1 0 0 2 2 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0  ]

Is it correct to say that for A we have 64 : 15 compression namely 4.2 : 1 cause we keep 15 coefficients out of 64? does the same apply for B ?


texts: filter coefficients to know high pass and low pass filter

texts:
Looked on few questions already there in dsp.stackexchange and couldn't get the direct answer. Basically I have a kernel of a gaussian separable filter.

float w[N*2+1] = {-0.00467444956f, -0.0531099625f, 0.0152470656f, 0.300147355f,  0.484779954f, 0.300147355f,  0.0152470656f, -0.0531099625f, -0.00467444956f};

float w2[N*2+1] = {-0.00899997167f, 0.0197822638f, -0.0300691985f, 0.0374629796f, 0.963647842f, 0.0374629796f, -0.0300691985f, 0.0197822638f, -0.00899997167f};

The coefficients in first first increasing and then starts decreasing.The middle point is strongest and it is close to 0.5.

In 2nd the coefficients are of alternate sign with middle being storngest close to 1.

What properties these coefficients signify?
How can I relate them to high pass or low pass filter?

I am using the above kernels in the smoothening the image using Gaussian separable filter.


texts: Connected Component Labeling - Equivalence Table

texts:
I am trying to write a connected component labelling program, I have done the first part to create a label based on the neighbours. For simplicity I am using a 4-connectivity. I am using ANSI C and to create the equivalence table I am using two linked lists with structures defined like this:

struct T2_node
{
    size_t value;
    struct T2_node *next;
};

struct T1_node
{
     struct T2_node *T2;
     struct T1_node *next;
};

Every newly found object has a T1_node which links to a T2_node linked list that keeps the labels of the objects that particular label is connected to. As an example you can see this thresholded image (very noisy!) I am trying to find the connected components of:

A section of the equivalence table for this image is:

label id: connected labels (including its self).

835: 847, 835,
834: 846, 834,
833: 820, 833,
832: 812, 832,
831: 831,
830: 830,
829: 818, 829,
828: 828,
827: 929, 778, 840, 783, 827,
826: 826,
825: 825,
824: 993, 944, 904, 916, 839, 824,
823: 823,
822: 957, 794, 837, 822,
821: 956, 912, 888, 848, 821,
820: 1021, 997, 915, 924, 887, 812, 833, 820,
819: 819,
818: 829, 818,
817: 817,
816: 816,

A grey scale image of the labels of the same image can be seen below:

As you see, the labelled regions within a connected component do not all touch!

My final question is: What algorithm should I use to effectively merge the labels? In other words, how can I find which labels belong to one connected object? I am a graduate student in Astronomy, so I do not know the details of graph theory and set theory or how to implement them in a computer program, take C for example! (although I am interested to learn). So if there is a specific method to do this, please also introduce a reference I can learn that method from, thanks in advance.


texts: What are the wavelet packet functions and how are they different from the wavelet and scaling functions?

texts:
I understand that in the decomposition process, wavelet and scaling functions are used to split the signal to approximation and detail coefficients.

Most of the time I saw the wavelet and scaling function associated with the normal discrete wavelet transform. I assumed the same functions are used for the discrete wavelet packets but now both the approximation and details are decomposed to form a binary tree.
The wavelet and scaling functions are set and have the same length regardless of the level of decomposition (based on the vanishing moment).
I was looking into wavelet packets since I want coarser decomposition and came across the wavelet packet functions

and now I am confused as these functions remind me of the continuous wavelet transform and the scaling process of the wavelet functions. They appear to be more scaled relative to the depth but I might be wrong and do not have the same length (based on the vanishing moments).
So what are these wavelet packet functions and how are they used for the decomposition process?


texts:
here is the Simulink model that i'm using :

The matlab function script :

function tau = fcn(sig,refsig)

tau = gccphat(sig,refsig);

The errors :


texts:
I am trying to prove if $h_1[n] = 2^{-n} u[n]$ and $h_2[n] = (0.1)^n$ are BIBO stable or not, but I'm not sure if these analysis are correct. I'm having some trouble to understand this specific subject.

For an impulse response $h[n]$ to be BIBO stable, equation $(1)$ need to be satisfied
$$
\sum\limits_{n=-\infty}^{\infty}\big\lvert h[n]\big\rvert = B <  \infty\tag{1}
$$

So, for $h_1[n]$,
$$
\sum\limits_{n=-\infty}^{\infty} \big\lvert 2^{-n} u[n]\big\rvert = \sum\limits_{n=-\infty}^{\infty} \big\lvert 2^{-n}\big\rvert\cdot\big\lvert u[n]\big\rvert = \sum\limits_{n=0}^{\infty} 2^{-n}
$$
according to Wolfram this sum is equal to $2$, and therefore, system $h_1[n]$ is BIBO stable.
For $h_2[n]$,
$$
\sum\limits_{n=-\infty}^{\infty} \big\lvert (0.1)^{n}\big\rvert =\sum\limits_{n=-\infty}^{\infty} (0.1)^{n}
$$
since $-\infty\le n\le\infty$, the system is unstable


texts: How does modulation frequency appear in the modulation spectra?

texts:
I'm trying to familiarise myself with the concept of modulation spectra that I met in this work. My main question appears in the title: how does the modulation frequency appear in the modulation spectra?
Calculation
It is clear that the modulation spectra is the power spectrum of the magnitude spectrogram of the STFT, while the word modulation spectrum is used to refer to the individual "spectrums" that the modulation spectra are composed of.
AM signal background
The expression for an AM signal is as a reminder:
$$ y(t) = (1 + \alpha \cdot\cos(2 \pi \omega_{m} t + \phi_{a})) \cdot \cos(2 \pi \omega_{c} t + \phi_{b}) $$
(changed due to suggestion below)
Motivation
However, it's not clear to me what's the motivation of this method. I think of the FFT as the frequency decomposition of a signal. Doing FFT twice would mean the inverse FFT for me (except for the scaling), but it is not the case because the phase is tossed away.
I can think of the spectrum as a time-domain signal, similarly as in the cepstrum approach where periodicities in the spectrum is of interest to find the harmonics, but that's also a different concept to this.
I don't understand where the modulation should appear and what are the quantities on the x and y axis of the modulation spectra.
My understanding so far
In the case of the AM signal, there should be three peaks, one corresponding to $\omega_{c}$, and the other to being $\omega_{m}$ distance away. In the previous version of my post, this didn't appear because of the 1 Hz frequency peaks were smeared by the mainlobe of the carrier frequency.
I also changed the way I calculate the modulation spectra. Each frequency bin can be thought of as a time-domain signal, so it makes more sense to take the FFT of the individual frequency bins signal, like here.
Example code
import numpy as np
import librosa
import matplotlib.pyplot as plt
import math
from scipy import signal

# - AM signal generation -
duration = 1 # in seconds
fs = 44100 # Hz
carrier_freq = 2000 # Hz
mod_freq = 800
x = np.linspace(0,duration,endpoint=True,num=fs)
audio = (1 + 1 * np.cos(2*math.pi*mod_freq*x)) * np.cos(2*math.pi*carrier_freq*x)

# - Visualise AM signal -
plt.plot(x,audio)
plt.xlabel("time (s)")
plt.ylabel("amplitude")
plt.xlim([0,0.02])

# - Calculate spectrogram   -
spectrogram = np.abs(librosa.stft(audio, n_fft=1024))

# - Spectrogram produced -
max_freq = fs // 2
plt.imshow(np.log10(spectrogram),aspect="auto",extent=[0,duration,max_freq,0])

# the carrier and the carrier - m/2 and the carrier + m/2 should appear
plt.xlabel("time (s)")
plt.ylabel("frequency (Hz)")
# - Calculate modulation spectra -
mod_spec = np.zeros((spectrogram.shape[0],spectrogram.shape[1]//2 + 1))

for bin in range(spectrogram.shape[0]):
    sajt = signal.windows.hann(spectrogram.shape[1])

    mod_spec[bin,:] = np.abs(np.fft.rfft(spectrogram[bin,:]))**2

# - Visualise modulation spectra -
# No idea about the axes!
plt.imshow(np.log10(mod_spec.T),aspect="auto",extent=[0,fs//2,fs//2,0])

Example figures

The plot seems to be correct.

From the spectrogram, I can identify the carrier frequency.

So, is the modulation frequency present here somewhere?


texts: Reducing sampling rate by a non-integer factor

texts: I have data that is sampled at 12 kHz. The downsample function in software such as MATLAB only allow you to downsample by an integer factor, i.e. from $12\to{6,4,3,2,1,\ldots}$ kHz. How would I downsample by a non-integer factor, say, from $12\to{10,9,8,7,5,...}$ kHz (pseudocode/description is fine)? You can assume that the signal has been appropriately lowpass filtered prior to downsampling. Also, the desired new sampling rate (after downsampling) is a "nice" number (i.e., I'm not trying to do $12 \to 9.7432$ or something like that).


texts: Filter response is changed according to its order

texts:
I want to build a band-reject filter for the 60 Hz noise. So, I tried to use Butterworth filter to implement this specification. However, when I chose an order of 8, I got the following response, which makes sense.

However when I chose an order of 30. I got this weird response. Why is that?

The code to get the coefficient for the 8 and 30 order respectively is shown below

# fs = 200
b, a = signal.butter(8, [58 / (fs / 2), 62 / (fs / 2)], btype='bandstop')
b, a = signal.butter(30, [58 / (fs / 2), 62 / (fs / 2)], btype='bandstop')


texts: Good features / algorithms for recognizing car-models in images

texts:
I have a question concerning object recognition, especially recognizing car-models! I am at the beginning of a work about identifying the same car-model in different images. At the moment I think one of the best algorithm for 3D object recognition is SIFT but after playing around a bit with a demo implementation I have the strange feeling this algorithm has some problems with shiny metal objects like cars, especially if they have different colors.

Does anyone know some work on this area in general some suitable algorithm for the task of finding the same car-model in different images?

Thanks in advance for your help!


texts:
I'm trying to plot 2D gaussian function using "imshow" in MATLAB. However I could not see any white points in image. I tried to change sampling values and limits of x and y values. In addition to this, I tried z.^2 and higher orders but nothing changed.Thank you for any help.What can be the problem? My code is below:

sigma = 10;
[X,Y] = meshgrid(-40:0.5:40,-40:0.5:40);
z     = (1/(2*pi*(sigma^2)))*exp(-((X.^2)+(Y.^2))/(2*(sigma^2)));
figure
surf(X,Y,z);%shading interp
figure
imshow(z);


texts:
As I found out in my previous question Where to get transtion matrix for Kalman filter?, I need a model for correct usage of Kalman filter.

In paper On the Intrinsic Relationship Between the Least Mean Square and Kalman Filters there is model-less Kalman filter - The Kalman filter for deterministic states:

\begin{align}
\textbf{g}(k) &= \textbf{P}(k)\textbf{x}(k) / (\textbf{x}^T(k)\textbf{P}(k)\textbf{x}(k) + r)\\
y(k) &= \textbf{x}^T(k)\textbf{w}(k)\\
e(k) &= d(k) - y(k)\\
\textbf{w}(k) &= \textbf{w}(k-1) + \textbf{g}(k)e(k)\\
\textbf{P}(k) &= \textbf{P}(k-1) - \textbf{g}(k)\textbf{x}(k)\textbf{P}(k-1)
\end{align}

where $\textbf{g}(k)$ is gain vector $\textbf{P}(k)$ is covariance matrix, $e(k)$ is error, $\textbf{w}(k)$ adaptive filter weights and $\textbf{x}(k)$ are inputs (not really a states here), $d(k)$ is target, $y(k)$ is filter output.

Nice thing is that $\textbf{x}^T(k)\textbf{P}(k)\textbf{x}(k)$ is scalar, so no matrix inversion is needed.

What is difference of this in comparison with RLS filter?
What should be better?
Is it case specific or is there some obvious difference what I am missing?
Is it even a good idea to use this Kalman filter over RLS or NLMS?


texts: how to understand the MIMO zero forcing as a bank of decorrelators?

texts:
I am studying the zero forcing as a simple linear receiver for MIMO detection. according to following notes:
http://www.ee.ic.ac.uk/bruno.clerckx/All.pdf.

The zero forcing can be understood as the following two ways.

For a single stream $\mathbf{q}$, we can project the received $\mathbf{y}$ onto a subspace orthogonal to channel matrix $\mathbf{H}$ without $\mathbf{h}_q$, then apply the matched filter to the projected signal or,
We can just work on all streams in one shot with the pseudo inverse of the channel matrix, i.e., $(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T$.

I can understand that the second approach certainly cancels the interference as $(\mathbf{H}^T\mathbf{H})^{-1}$ exists if columns of $\mathbf{H}$ are independent, and obviously $(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{H} = \mathbf{I}$ such that $(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{y}=\mathbf{x} + (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{n}$.

It seems to me that according to (2), a decorrelator for the $q^{th}$ stream is just the $q^{th}$ row of the pseudo-inverse of the matrix $\mathbf{H}$.

My question is, can we prove, or conclude (2) from (1)?

The following is the MATLAB script implemented according to the methods described in the above-mentioned notes.

%zero forcing

randn('seed', 0);
m=2; n=2;

x = randn(n,1)+1j*randn(n,1);
nz=0;

H = randn(m,n)+1j*randn(m,n);

y=H*x+nz;

%zeero forcing filter

G0 = inv(H'*H)*H';
x0=G0*y;

%disp(x-x0);

% transform the y into the null space of the

% subspace spanned by H1m (H w/o h1)

h1=H(:,1);
H1m=H(:,2:end);
[U, S, Vc]=svd(H1m);

% Q1 spans null space of the subspace spanned by H1m (H w/o h1)

% (actually it is null space of H1m')

% i.e., Q1'*y projects y into the null space of H1m

% y1 = Q1'*y = Q1'*h1*x1 + Q1'*n + zero (cancelled interference)

Q1=U(:,n:m);

% (Q1'*h1)' is the match filter to the transformed y1

% g1=h1'*(Q1*Q1'), note that Q1*Q1' is not orthogonal but Q1'*Q1=I

g1=(Q1'*h1)'*Q1';

% normalization

g1=g1/(norm(g1)^2);
x01=g1*y;

%x(1)-x01

%snr for stream 1

snr1=(norm(Q1'*h1))^2

nv=inv(H'*H);

%snr1 - 1/nv(1,1)
```


texts: Demodulating upper sideband (USB) signals

texts:
A modulating signal $m(t)= \sin(2 \pi \nu_m t)$ is transmitted via a carrier of frequency $\nu_c$ using upper sideband (USB) modulation. Show that the USB modulated signal can be demodulated using a reconstructed carrier signal and a low pass filter.

Attempt:

USB (SSB signals in general) can be demodulated by multiplying them with a reconstructed carrier $\cos(2\pi \nu_c t).$ So I first tried to write an expression for the transmitted USB modulated signal:

$$f(t)=A \Big[ m(t) \cos(2 \pi \nu_c t) - \hat{m}(t) \sin(2\pi \nu_c t) \Big]$$

$$f(t)=A \Big[ \sin(2 \pi \nu_m t) \cos(2 \pi \nu_c t) - \hat{m}(t) \sin(2\pi \nu_c t) \Big] \tag{1}$$

Where $A$ is some amplitude, and $\nu_c >> \nu_m$. So to proceed further, I think we need to calculate the Hilbert transform $\hat{m}(t).$ It is possible to do this by finding the imaginary part of the analytic signal:

$$m_a (t) = m(t) + j \hat{m} (t)= 2 \int^\infty_0 M(\nu) \exp(j 2 \pi \nu t) d\nu \tag{2}$$

But the Fourier transform of $m(t)$ is tedious, so I am wondering if there is a simpler way of doing this?


texts: Direct downconversion (IQ downconverter) with ASK modulation-- what kind of Balun do you need?

texts:
I have a simple but confused question.
I am trying demodulate a 915MHz ASK modulated signal with direct downconversion. So what I am building up is:
Antenna --> SAW bandpass filter @915MHz --> LNA --> Direct IQ down-converter (max2021 datasheet LO = 915MHZ) --> balun? --> decoder(sdr)
The max2021 downconverter is a 4-ports output: I+,I-,Q+,Q-. To put the differential signal back to single-ended, I was taught this is a balun's job.
But when I trying to build this system the balun concept is confusing. If successfully down-converted and DC blocked  , the IQ outputs are just unmodulated differential square pulses-- so what kind of balun do I use and at what frequency? Do I use a digital buffer in this case instead of an analog balun?
ps: I am using direct downconversion for carrier interference cancellation(move the carrier interference to DC then DC block) and design simplicity; and IQ downconverter to avoid phase ambiguity.


texts: Moving average filter magnitude attenuation

texts:
I am dealing with a moving average filter which filters power system signals to rid them of noise or harmonics etc. The fundamental or useful component of these signals is centered around 50 Hz or 60 Hz. The magnitude response of the filter at an angular frequency say  as per some standard papers is given as
$|H(e^{(j\omega_1)})|=|\frac{\sin(L_m \omega_1 /2)}{L_m \sin(\omega_1/2)}|$ ,
 where $L_m$ is the window length of the filter and the value of $\omega_1=\frac{2\pi f_f}{F_s L_m}$ . Sampling frequency is given by $F_s$
whereas the input signal frequency to the filter is given by $f_f$ .  The idea is when the filter input signal frequency is equal to the fundamental (no noise or harmonics etc)  or even close to it (say 59.5 Hz or 59.7 Hz considering a 60 Hz system) the magnitude response (gain) of the above mentioned filter should be equal to one  or close to it for any sensibly chosen value of the window length and the sampling frequency. For example if i choose a value of $L_m$=31 and and consider the input signal frequency $f_f$ to be 59 Hz, with a sampling frequency $F_s$  of say 1440 Hz, i end of up with a gain of 0.9972, then if i use $f_f$=58, i end up with a gain of 0.9973, similarly along the lines, for $f_f$  equal to 50 Hz i get a gain of 0.9980, showing that gain increases as frequency deviates from the nominal (which does not make any sense). I am having a hard time understanding this, could you clear this confusion for me. Could you show me how the gain is close to unity for small frequency deviations from the nominal


texts: Why Low Frequency called Approximation and High Frequency Detail?

texts: While we are decomposing a signal using Wavelets into levels why do we call our high frequency components as Detail and Low Frequency Components as Approximation ?


texts: Cascading first order digital filters in C++

texts:
In a related question a probable solution was given to build a first-order digital filter and then cascade three of them in order to turn white noise into pink.  I have applied the C++ as follows but still the signal sounds and looks like white noise.
I would like to know what is wrong with my implementation and why I don't hear pink noise.  For reference, the poles and zeros come from Robert Bristow-Johnson's work here.
Header:
float *state = nullptr;
Implementation:
state = new float[0.0]; in Constructor.
Then in the loop, for (int i=0; i < numSamples; i++)
float first = first_order_filter(whiteNoise, 0.99572754, 0.98443604, state);
float second = first_order_filter(first, 0.94790649, 0.83392334, state);
float third = first_order_filter(second, 0.53567505, 0.07568359, state);
out1 = third;

Where first_order_filter is defined as in Robert's answer here:
https://dsp.stackexchange.com/a/70963/11391
I would love to know if this code is approximately correct/ where the problem lies.


texts: Unable to understand SNR equation in Matlab

texts:
Can someone please help me understand this equation of SNR given in MATLAB:

SNR_dB = 1:0.1:7
SNR = 10.^(SNR_dB/10)


texts: Python: Least Squares Support Vector Machine (LS-SVM)

texts:
I'm looking for a Python package for a LS-SVM or a way to tune a normal SVM from scikit-learn to a Least-Squares Support Vector Machine for a classification problem.

The goal of a SVM is to maximize the margin while softly penalizing points that lie on the wrong side of the margin boundary. The function that is used is a Quadratic Programming (OP) problem.

A LS-SVM which defines a least squares cost function and replaces the inequality contraints with equality constraints and is a Linear Programming (LP) problem.

Wikipedia: LS-VSM

Does someone ever tried something like that and can help me?


texts: Difference between a digital lock-in amplifier and a FFT when extracting phase of a signal?

texts:
I am trying to measure the relative phase of a sine wave fixed at a particular frequency in a noisy environment. My initial approach is to simply collect $N$ samples, take an FFT, and then just extract the phase at the operating frequency (which I know apriori).

When reading about methods to extract signals from noisy environments, I came across the lock-in amplifier. However, I am confused about whether I should have any reason to expect it to perform better than the FFT (ignoring run time, just looking at ability to extract my phase).

Specifically, imagine that I implement the lock-in amplifier digitally. I would then numerically calculate the following two integrals (source):

\begin{align}
X_{LI} &= \frac{1}{T}\int\limits^{T}_{0} U_{\rm in}(s) \cos(2\pi f s) \,ds\\
Y_{LI} &= \frac{1}{T}\int\limits^{T}_{0} U_{\rm in}(s) \sin(2\pi f s) \,ds
\end{align}

and then I would get my phase using:

$$
\theta_{LI} = \tan^{-1}\left(\frac{Y_{LI}}{X_{LI}}\right)
$$

However, if I were to use an FFT, the FFT would evaluate

\begin{align}
X_{FFT}(f) + iY_{FFT}(f) &= \int\limits_0^T U_{\rm in}(s)e^{-2\pi i s f}\,ds\\
&=\int\limits_0^T U_{\rm in}(s)\cos(2\pi f s)\,ds + i \int\limits_0^T U_{\rm in}(s)\sin(2\pi f s)\,ds
\end{align}

Which seems to imply that just as with the lock-in amplifier, I would get

\begin{align}
X_{FFT} &= \int^{T}_{0} U_{\rm in}(s) \cos(2\pi f s)\,ds = \frac{1}{T}X_{LI}\\
Y_{FFT} &= \int^{T}_{0} U_{\rm in}(s) \sin(2\pi f s)\,ds = \frac{1}{T}Y_{LI}\\
\theta_{FFT} &= \tan^{-1}\left(\frac{Y_{FFT}}{X_{FFT}}\right) = \tan^{-1}\left(\frac{Y_{LI}}{X_{LI}}\right) = \theta_{LI}
\end{align}

So does this mean that I would always get the same result whether I use a digital lock-in amplifier, or a digital FFT  on the same data set?
What is the benefit of using a lock-in amplifier?
Is there no advantage to the lock-in amp in my application?
What is then an example of an application that the lock-in amplifier is best suited for?


texts: Interpreting FFT Phase - why phase of a cosine?

texts: I understand how to interpret the magnitude result from the FFT, but why is the phase that we obtain, arctan(Im(x)/Re(x)) indicative of the phase shift of a cosine graph, and not a sine graph?


texts: What is the power spectral density after filtering?

texts:
Consider a signal $x(t)$, which is input to a pulse shaping filter with transfer function $g_t(t)$:

$$x(t) = \sum_n d_n  \delta(t-nT_s)$$

with $n$ an index from negative to positive infinity and $d_n$ random, equiprobable binary symbols such as 1 and -1 (antipodal) or 0 and 1 (on-off keying).

Well, the power spectral density of the filtered output is given by:
$$S_{ss} = \frac{1}{T_s}  E\left[|d_n|^2 \right] |G_t(f)|^2 $$
where $E[\cdot]$ denotes the expected value.

I'm confused by this formula since up to now I thought the power spectral density is given by the signal mutliplied with the absolute value squared of the transfer function. But why do I have to consider the expectation value in this case and why am I scaling by the sampling time?

Thx for any help!


texts: Additive Gaussian white noise - bandwidth

texts:
The average power at output when additive Gaussian white noise is present is given by the formula below:

$$P_0 = \int\limits_{-B}^{B} \frac N2\ df ,\quad \text{with  $B$ being the bandwidth.}$$

My question is when the bandwidth is in terms of $\textrm{rad/s}$, what would be the formula?
The one stated divided by $2\pi$ or another formula?


texts:
We are looking at FMCW radar ADC signals and converting from the time domain to the frequency domain in order build a range doppler map but want to have the best signal to noise ratio possible.

In order to reduce spectral leakage, we are using the Hann function.  It just seems to be accepted that we use the Hann function.

Are there any other functions that we should be considering?


texts: Applying multiple IIR and FIR filters at once

texts:
I would like to appy multiple filters to a signal at once.
Let's assume I got two frequency bands I want to extract. I compute a Butterworth Filter for 300-500Hz and one for 700-800Hz. I would like to combine those filters into a single filter and apply it to the signal. The resulting signal should contain only frequencies in range 300-500Hz and 700-800Hz.

I assume that if the filters are FIR (finite impulse response) filters I just need to add them. Using that single filter, I apply a convolution to the signal. Is that correct?

I am pretty lost if it gets to IIR (Infinite impulse response) filters which are - as I understand - divisions of polynoms in a radial space I do not understand. How can I combine two IIR filters?

Kind regards,
Thomy800