texts: Removal of noise from vibration recordings

texts:
Problem
I am analyzing accelerometer recordings that are contaminated by 'spikes' caused by issues in the digital circuitry. The spikes are problematic, because I wish to determine the peak-amplitude from these signals using an automated procedure from a large number (hundreds, perhaps thousands) of these recordings. Hence, manually determining the amplitudes is possible, but not preferred.

Background
The acceleration recordings are oscillatory in nature and reach frequencies of ~200 Hz. I am planning on using a low-pass filter in Matlab to remove the contaminating spikes in the signal. Sample rate is around 730 Hz. An FFT-based filter may be complicated by the fact that the signal gradually increases and decreases in frequency over time.

Here is a picture of an example recording including the pesky spikes (the noise) in the first 100 ms:

Question
I am using Matlab for signal processing and I am looking for advice on which digital filter for use as a high-pass filter? For example, Butterworth and Chebyshev filters are quite commonly applied if I am not mistaken, but why? Will wavelets do me any good?


texts: Has there been any research on time-windowed spectral estimation using more than one window function at the same time?

texts:
In the Welch's method, a signal is divided into overlapping segments, multiplied with a window function, and transformed to frequency domain via fourier transform.

Why does the method use only one window function? Suppose I would use (say) ten window functions to produce ten spectral estimates. By averaging the estimates, my SNR would be higher than by using just one of the windows?


texts: Mapping ring asteroids in Elite Dangerous by image processing

texts:
In the game Elite Dangerous, which is an online spaceflight-based sandbox game set in a procedurally-generated galaxy, one of the most profitable in-game occupations is mining commodities from asteroids in the rings of certain planets.  Because mining can be so profitable, it is widely used to gather in-game currency to fund other activities.  Many players do not actually enjoy mining, but engage in it solely to fund their preferred gameplay, so Credits per Hour are very important to them.

The content of the ring asteroids is constant,  as is their position in 3-dimensional space over time, however many rotate in-place at a constant rate around a constant axis.  Rings are effectively circles with no centre thousands of kilometers wide, and are around ten kilometers thick.

I believe that the pattern of asteroids in the rings may repeat, though the mineable content of the asteroids may not repeat, or may repeat over multiple repeats of the basic 'tile(s)'... Determining that is part of the point of this project.

Elite Dangerous logs many things that are encountered in-game, but ring mining activities are not logged, nor is there any coordinate system accessible to the players.

The mining process involves dropping out of supercruise (a rapid travel mechanism) near a ring, flying to that ring, and identifying suitable asteroids to mine, either by randomly prospecting nearby asteroids until suitable asteroids are found, or by pre-mapping high-yield asteroids by the expedient of screenshots, and using them in conjunction with fixed points of reference to orient oneself and find and mine each successive high-yield asteroid. There are a limited number of asteroid shapes, and all asteroids in a ring are similarly coloured, except when a pulse-wave scanner is used, which turns certain asteroids different colours according to their content.

Once mapped, a player could drop into a mapped area, and the companion software would highlight each asteroid according to its status, and point the player to the next most profitable asteroid.

So, that being the background, this is the problem that I wish to solve: I would like to write a companion application to Elite Dangerous that uses image processing techniques and optical character recognition to screen-scrape the relevant ring metadata.  Whenever a player running this companion app drops into a ring, it should use photogrammetry to locate the player within an existing map or begin constructing a new map.  Whenever an asteroid is prospected, the app should be able to log the asteroid's content and associate it with the relevant asteroid.

In order to write this application, (using C#.NET), I need to know the best approach to using photogrammetry to map a ring in as close to realtime as possible.  Are there any algorithms or APIs that would be particularly suitable for this task?  Do any of them lend themselves to GPU processing?  I would prefer free open source solutions, given that this will most likely be freeware.

I don't need to discuss OCR... just photogrammetry of the ring asteroids.

Here is a video that includes footage of the rings in question from about 3' 35".

Here's a Google Drive folder-full of mining guide photos.

For example:

The algorithm need not map the shape of each asteroid in any great detail, it need only place each asteroid reasonably accurately within 3d space.  The shape and position of asteroids is constant,with the exception of rotation, and the content of asteroids can be identified by position alone.  As some asteroids rotate very slowly, and are highly irregular in shape, attempting to identify the shape of the asteroids would be unnecessary work.

Additionally, since edge detection was mentioned, Elite Dangerous has a "night vision" mode which functionally highlights the edges of objects.  It is entirely acceptable to require players using the proposed application to activate that feature, since much mining is done with it activated anyway.


texts: Discrete implementation of the PI controller

texts:
I have been implementing discrete PI controller in the incremental (velocity) form in C++.
I have been looking for the anti-windup mechanism. One idea which I have can be described by
following pseudocode:
    integral_gain = proportional_gain*sampling_period/(2*integral_time_constant);

    // e(k)
    error = reference_value - actual_value;

    // dup(k)= kp*[e(k) - e(k-1)]
    proportional_increment = proportional_gain * (error - error_previous);
    // trapezoidal integration rule i.e. dui(k) = ki*[e(k) + e(k-1)]
    integral_increment = integral_gain * (error + error_previous);
    // du(k) = dup(k) + dui(k)
    action_increment = proportional_increment + integral_increment;
    // u(k) = u(k-1) + du(k)
    tmp = action_previous + action_increment;

    // anti-windup
    if (tmp > action_max) {
      action = action_max;
    } else if (tmp < action_min) {
      action = action_min;
    } else {
      action = tmp;
    }

    // u(k-1) = u(k)
    action_previous = action;
    // e(k-1) = e(k)
    error_previous = error;

In math:
Let $x_k$ be the "action".  Further, define an anti-windup function as
$$f_{aw}(x) = \begin{cases}
x_{max} & x > x_{max} \\
x_{min} & x < x_{min} \\
x & \mathrm {otherwise}
 \end{cases} $$
Then calculate $x_k$ from:
$$
\Delta p = k_p (e_k - e_{k-1}) \\
\Delta i = k_i (e_k + e_{k-1}) \\
x_k = f_{aw} \left ( x_{k-1} + \Delta p + \Delta i \right)
$$
I am not really sure whether the anti-windup mechanism mentioned above will be appropriate despite the fact that the incremental form is used. Can anybody help?


texts: Can anybody tell me why do we use single pole filter? And may be an example of Single pole filter in matlab. Good literature recommendations are also welcomed.Moreover, that would be good if anyone graphically explains, what happens to the signal when single pole filter is applied?


texts: Using IFFT for the purpose of "Multiple Orthogonal Frequency Shift Keying"

texts:
Is it possible to create a set of complex orthogonal sinusoids using FFT or IFFT
In other words a kind of  "Multiple Frequency Shift Keying" using  multiple sinusoids that are created with a fourier transform ( FFT or IFFT) and where one frequency is sent at a time?


texts:
I was truing to solve an example of DTFT which is following multiplication property. The problem is
$$ a^n \sin(\omega_0 n) u[n]$$ we know that the definition of DTFT is
$$ X(j \omega) = \sum _ {n=-\infty} ^ {{+\infty}} x[n]e^{-j \omega n}$$
Multiplication in Time domain will be convolution in DTFT.

If we take the DTFT of $a^n u[n]$ we have $$\frac {1}{1-ae^{-j \omega}}$$ and DTFT of $\sin(\omega_0 n)\,u[n]$ will be  $$\frac{\pi}{j}\sum _ {l=-\infty} ^ {{+\infty}}  \delta(\omega + \omega_0 - 2\pi l) - \delta(\omega - \omega_0 - 2\pi l)$$

I have confusion how can I write it in the form of multiplication property.


texts: Applying Image Filtering (Circular Convolution) in Frequency Domain

texts:
To filter an image we can:

Use a 3x3, 5x5, 7x7, etc. filter, that is convolve the image and the filter in the space domain.
Use a FFT on both the image and the filter, multiply them together in the frequency domain, and then perform an IFFT on the image to get the filtered result.

These two approaches should produce equivalent result.

Questions:

Are there any conditions placed on the filter so that the results are equal?
Will there be any information loss in the second approach or will the images be exactly equal (neglecting floating point precision error)?
As the kernel size increases, will the second approach start outperforming the first?


texts: PIsarenko Harmonic Decomposition method for PSD estimation

texts:
i would like to present  steps for estimation Power Spectrum density using Pisarneko Harmonic  Decomposition

i would like to present code for  step 1 part

function [vmin,sigma]=phd(x,p)

% the Pisarenko Harmonic decomposition
x=x(:);
R=Covar(x,p+1);
[v,d]=eig(R);
sigma=min(diag(d));
index=find(diag(d)==sigma);
vmin=v(:,index);
end

now let us  consider following time series

fs=100;
ts=1/fs;
t=0:ts:2.93;
x=10*sin(2*pi*10*t)+15*cos(2*pi*10.5*t)+14*sin(2*pi*10.19*t)+3*randn(size(t));

PSD using music method can be estimated  as follow

pmusic(x,6);

how can i  got the same result using pisrenko harmonic decomposition? thanks in advance


texts: Obtain the frequency response of a time domain chirp signal

texts:
I have simulation data in the time domain and I need to represent them in the frequency domain.

The data have been obtained time integrating the following equation:

$m \ddot{x} + c\dot{x} + kx = A \sin{(\phi(t))}$

where $\sin{(\phi(t))}$ represents a logarithmic chirp. Therefore, since also the integrated time histories are chirps, I cannot use the FFT to transform them in the frequency-domain, but I have to do something like it was done in this answer: system response: time vs frequency. Why do I get different magnitudes?

Therefore I used the spectrogram in this way:

% Resample at a constant sampling frequency
dt = 0.0001;
t0 = time(1);
tend = time(end);
new_t = t0:dt:tend;
Fs = 1/dt;

acce_TS = timeseries(acce_cog,time);
acce_res = resample(acce_TS,new_t);

% Spectrogram
win = hamming(256);
noverlap = 0;
nfft = 512;
X = abs(spectrogram(acce_res.data, win, noverlap, nfft,Fs));
X = 2*X/sum(win);
surf(X)

And here come the problems:

1) On the time axis I have wrong values: the time history lasts $49.5 \ s$ whereas in the spectrogram plot I see values from 0 to 2000

2) The peak response should be at nearly $30 Hz$ (which is the resonance frequency of the system, given as $f_n = \sqrt{k/m}$), whereas from the spectrogarm I see it at nearly $5 Hz$

3) The amplitude of the peak of the spectrogram is a bit higher than in the time history

4) I actually need to see only the frequency-amplitude plot, not the 3D representation, so how can I do that? At the moment I am using view([90 0]), but is there a way to plot directly the 2D plot?

EDIT: Here the plots


texts: Removing transients in highpass filtering with MATLAB

texts:
I've made a simple first order IIR highpass filter with a zero at z = 1 and a pole at z = 0.9. Its frequency response looks like this:

Now, I filter a DC signal using this filter. Here's the MATLAB code I use to do it:

b = [1 -1]; % Zero at z = 1
a = [1 -0.9]; %Pole at z = 0.9

figure(1)
freqz(b, a)

t = 1:100;
x(1:length(t)) = 1; % Constant function

y = filter(b, a, x);

figure(2)
plot(t, x)
xlabel('Time');
ylabel('Input Signal');

figure(3)
plot(t, y)
xlabel('Time');
ylabel('Output Signal');

As my filter is highpass, I expect the DC to become zero, or atleast become severely attenuated. However, the output signal I get looks like this:

From my understanding, this exponential output is a transient produced because I haven't set the initial conditions correctly. Sure enough, setting x[-1] = 1 solves the problem. However, this works only for this particular input DC signal. For any general input signal, how do I set the initial conditions so that transients aren't produced?

Edit : I'm aware that the filtfilt() function does forward-backward filtering with transient minimization, but I really want to port the filter to an embedded platform, so I need to understand how transient removal works. Thanks in advance for the help!

Edit 2 : As suggested by Kuba Ober, I tried setting x[-1] as the value that it actually should have been. It works fine for a DC input, but here's what happened for a sinusoidal input:

clc; clear all;
p = 0.9

a = [1 -p]
b = [1 -1]

n = 1:100; % Samples
f = 0.2; % Frequency in Hz
Fs = 10; % Sampling rate in samples per second
t = n/Fs; % Time axis

x = sin(2*pi*f*t);

% Filter with the appropriate initial conditions
y = filter(b, a, x, filtic(b, a, [], [sin(2*pi*f*0)]));

figure(1)
plot(t, x)
xlabel('Time');
ylabel('Input Signal');

figure(2)
plot(t, y)
xlabel('Time');
ylabel('Output Signal');

Here's the input signal :

And here's the output :

The first peak is visibly smaller than the second, which indicates some transients being present. I'm not entirely sure about this, but I think the reason it doesn't work is because just setting x[-1] is not enough, I also need to set y[-1]. The problem here is that there's no way to find out what y[-1] actually should be.

Edit 3 : Let me provide a little more info on the problem I'm working on. I'm trying to use filters to remove noise from ECG (Electrocardiogram) signals in an embedded platform. Here's a typical ECG signal, after filtering:

Here's what an ECG signal looks before filtering:

Note the DC offset in the signal before filtering. For filtering, I need a notch filter to remove high frequency power line noise and a highpass filter to remove the DC and the low frequency "drifting" of the signal.

The filters I use need to be linear phase, since the time domain morphology of an ECG signal is very important for diagnosis. However, my filter doesn't need to be real-time, as I'm doing the processing offline after acquiring the ECG signal from the patient. So, for implementing nonlinear phase IIR filters, I'm currently using forward-backward zero phase filtering.

One opinion that's shared by @Matt L. and @Drazick is that transients are unavoidable in real-time filtering and that I should use a longer input signal and crop off the first few seconds of the filtered output instead. This is something I'd like to avoid, as acquiring long ECGs from a living patient is somewhat difficult. Also, I do not have to filter in real-time, so any technique of transient removal that hinges on knowing the entire signal in advance is perfectly all right. Any help is appreciated!


texts: Power Spectral Density of Brownian Motion despite non-stationary

texts:
Note: I originally asked this on Physics Stack Exchange but haven't attracted any interest there so I'm asking here where it may be more relevant.

A white noise process, $\xi(t)$ with delta correlated two-correlation function $\langle \xi(t_1)\xi(t_2)\rangle = \delta(t_1-t_2)$, is clearly stationary  and has a power spectral density which is the Fourier transform of the auto correlation function a la Wiener-Khintchine theorem.

$$
S_{\xi \xi}(f) = \int e^{i 2 \pi f t} \langle \xi(t)\xi(0) \rangle dt = 1
$$

We call such a process white noise because it has a flat power spectral density as indicated by the above equation.

I am concerned here with the power spectral density of a Brownian motion process or a Wiener process, $W(t)$ which is the integral (in some sense) of a white noise process.

$$
W(t) = \int_{t'=0}^t \xi(t) dt = \int_{t'=0}^t dW(t')
$$

See the Wikipedia article on Brownian noise. In that article it is pointed out that since Brownian motion is the integral of white noise we have (using the Fourier transform convention above)

$$
\mathcal{FT}[W](f) = \frac{1}{-i 2\pi f} \mathcal{FT}[\xi](t)
$$

Then the argument follows that the power spectral density is like the absolute value of the Fourier transform so we should have something like

$$
S_{WW}(f) = \frac{S_{\xi\xi}(f)}{|-i 2\pi f|^2} = \frac{1}{(2\pi)^2} \frac{1}{f^2}
$$

That is, the Brownian process has a $\frac{1}{f^2}$ spectrum, or Red or Brown noise.

We can also perform a direct computation of the power spectral density from the definition of power spectral density as the limit of the magnitude of the windowed Fourier transform.

$$
S_{WW}(f) = \lim_{\tau \to \infty} \frac{1}{\tau} \int_{t_1=0}^{\tau} \int_{t_2=0}^{\tau} e^{i 2\pi f (t_1-t_2)} \langle W(t_1)W(t_2) \rangle dt_1 dt_2
$$

It can be shown for a Wiener process that the two-time correlation function is $\langle W(t_1) W(t_2) \rangle = \text{min}(t_1,t_2)$. We can write the integral as

$$
\frac{1}{\tau}\int_{t_1=0}^{\tau}\int_{t_2=0}^{t_1} e^{i 2\pi (t_1-t_2)} t_2 dt_1dt_2 + \frac{1}{\tau}\int_{t_1=0}^{\tau}\int_{t_2=t_1}^{\tau} e^{i 2\pi (t_1-t_2)} t_1 dt_1dt_2
$$

Mathematica (or a diligent calculus student) evaluates this integral to

$$
\frac{2\pi f\tau - \sin(2\pi f \tau)}{4 \pi^3 f^3 \tau}
$$

In the limit of $\tau \to \infty$ the second term vanishes and we are left with

$$
S_{WW}(f) = \frac{1}{2 \pi^2} \frac{1}{f^2}
$$

This is missing a factor of 2 compared to the formula above. I think the reason is related to the fact that I've only considered positive time for the windowed Fourier transform in the second derivation.

Here is my question
I have given two sketches that show Brownian motion has power spectral density $\sim \frac{1}{f^2}$. However, I have also stated the two-time correlation function for Brownian motion is $\langle W(t_1)W(t_2) \rangle = \text{min}(t_1,t_2)$. This shows that Brownian motion is clearly not a stationary process. I know that the Wiener-Khintchine theorem tells us that if a process is stationary then the power spectral density is given by the Fourier transform of the two-time correlation function. But I also had the impression that non-stationary processes didn't have well-defined power spectral densities, or that you could only define a time-windowed power spectral density. However, Brownian motion seems to be an exception that rule. It seems to be an example of a process which is non-stationary (so you can't use the Wiener-Khintchine theorem, for example) yet still has a well-defined power spectral density. Am I correct in this characterization? Can anyone give me more information about non-stationary processes that still have power spectral densities?


texts:
I was just wondering if there were any "obvious" resources, papers, or textbooks I should look into regarding the nyquist-criterion in higher dimensions. For example, this would be relevant if there is some curve that is i.e. not 1-to-1 w.r.t. time: $f: t \in \mathbb{R}^1 \rightarrow f(t) \in \mathbb{R}^1 $, but instead for some other parameterization, we have a line $l: u \in \mathbb{R}^1 \rightarrow l(u) = \big(\tilde t(u),\tilde f(u)\big) \in \mathbb{R}^2$, it is 1-to-1, but not onto, as $(\tilde t,\tilde f)$ is on a 1-dimensional manifold (a line) in a 2-D plane.

Then I'd like to know what sampling scheme I need to be able to recreate the curve through reconstruction. Using this approach, you could use nyquist's theorem for general spatial surfaces and manifolds as well.

This might be a math question, so I can migrate this. But I figured some experts on this sub-stack might know.

EDIT: several small edits for clarity/correctness.

EDIT 2: Fixed parameterization!!!

EDIT 3: I want to emphasize that although I am parameterizing by $t$, it may be more useful to think of $t$ as a spatial dimension, since time can only move forward.


texts: PLL carrier recovery: intermediate versus baseband

texts:
When looking at carrier recovery, a PLL is often applied to correct the frequency offset between the local oscillator and the transmitter after down conversion.

I'm wondering what happens to the PLL if this frequency offset is small? Does the PLL break as the carrier tone isn't as obvious and the loop filter can't as easily differentiate the modulation noise from the carrier?

Wouldn't it be more advantageous to have a local oscillator down-convert the signal to an intermediate frequency and then have the PLL track that?


texts: Connections between quantum field theory and signal processing

texts:
I've been working with quantum field theory on the basis of it being a stochastic signal processing formalism for a number of years. The projection to fourier components within the forward light-cone in 1+3 dimensions is comparable to projection to positive frequency components in 1+0 dimensions (and hence comparable to projection to the analytic signal).

Modern physics is largely about storing and processing digital records of multiple (often multiple thousands of) electronic signals. The digital records are most often somewhat ad-hoc, very sparse, and very nonlinear in their relation to the original signals, typically as times at which an analogue voltage makes a transition between distinct metastable levels ("a particle was detected at time ... and position ...", but in fact, for example, a micron-sized device transitioned to a current-flowing state for a few nanoseconds before being reset to a zero current state), nonetheless this is signal processing.

As a partial motivation for thinking this way, we can construe quantum fluctuations, in a signal processing PoV, as a Poincaré invariant noise, very comparably to the translation and euclidean invariance of thermal noise (see my Phys. Lett. A 338, 8-12(2005); non-paywalled arXiv preprint).

I'm very much aware of Leon Cohen's work on relationships between deterministic SP and quantum theory (but very happy to hear of anything recent in the SP literature), and Hilbert spaces are routinely used for deterministic SP, but I'm interested to hear of any uses, in the literature or otherwise, of Hilbert space mathematics for stochastic SP (particularly any that take an SP approach, but stochastic fields are also called "random fields" in the mathematical literature).

I'm also interested in Hilbert space methods for nonlinear control theory that take a more-or-less SP approach, if anyone can point me to such (particularly stochastic nonlinear control theory, but deterministic if needs must).


texts: Permutation test on two simultaneously recorded EEG signals

texts:
I have two EEG signals that are simultaneously recorded (and has been further synced) using two EEG devices, namely device 1 (ground truth) and device 2 (investigational device).
Imagine one device recorded from channel F4 and the other from Fp2. Now my aim is to validate that the investigational device works as accurate as of the ground truth and to do this I used different measures, namely: PSD and window-by-window Pearson correlation (over all subjects).
To generalize my results, I wanna apply a permutation test. But I am not sure about the null and alternative hypotheses I should define. Below I show the ones that I currently have (which is usual definition of H0 during permutation test):
H0: No difference between investigational and ground truth EEG recorders.
H1: They differ
But I wanna check whether the dependency between signals from device 1 and 2 are based on chance or the difference is significant. So, to conclude that if I find a p-value<0.05--> the similarities are not due to chance and if the p-value > 0.05 --> recording devices differ.
Can you please help me with redefining the H0, H1 e.g. for comparing PSDs?


texts:
Is it possible to realize DFT of 6 samples by using 3 DFT2 blocks? For example if had 2 * N samples by one step of decimation in frequency or time is possible to reduce DFT(2N) to 2 * DFT(N). So if we had 2 * DFT(N) blocks we could do it that way. What if we had 3N samples and 3 DFT(N) blocks. In my question N equals 2. How would the structure look and what is formula for X(k) element of spectrum?

The basic problem here is how to find different ways to calculate DFT with given blocks and using only additional sum blocks and multiplying with twiddle factors?


texts: Same z transformed function, but different answers of inverse z transform?

texts:
Given a $\mathcal{Z}$ transformed function $E(z)=\frac{1}{z+4}$.

I know there are several ways to get the inverse $\mathcal{Z}$ transform of this function :

Using partial fraction

$$E(z)=\frac{1}{z+4}$$

$$\frac{E(z)}{z}=\frac{1}{z(z+4)}$$

apply partial fraction here,

$$\frac{E(z)}{z}= \tfrac14 \left(\frac{1}{z} - \frac{1}{z+4} \right)$$

so $E(z)$ is,

$$E(z)= \tfrac14 \left(1 - \frac{z}{z+4} \right)$$

as you know, it is easy to use inverse Z-transform here.

$$ e[n] = \tfrac14 ( \delta[n] - (-4)^n ) $$

inversion formula method

$E(z)$ has simple pole at $z=-4$, the residue is evaluated as

$$\begin{split}
\text{residue}_{z=-4} &= ((z+4)\times E(z)\times z^{k-1})_{z=-4}\\
&= (1\times z^{k-1})_{z=-4}\\ &= (-4)^{k-1}\\ &= -(1/4)\times (-4)^k \end{split}$$

The problem is that inverse $\mathcal{Z}$ transform of the same function $E(z)$has 2 different answers. As you see, there is no delta function when I used inversion formula method.

At first, I thought it was because a simple pole, $z=-4$, is outside of the unit circle, but same things happen when I used a simple pole which is inside of the unit circle.

Why does this difference happen? What is the real answer?


texts:
I am doing digital signal processing for a software defined radio. I have a signal coming in packets of data (several samples) at an average sample rate of approximately some known value FSin. I am given a request for packets to be sent out periodically, with an average rate out of approximately FSout. If I knew these rates exactly I could resample my data by a ratio R=FSout/FSin and use a FIFO to make sure data always has somewhere to go when I get a packet and there's always data to send when I'm asked to send a packet. However, the clocks of the radio and my computer are not perfectly synchronized, so there is no exact known ratio for R. What would your solution be?

This is what I am thinking: I will have to adjust the ratio by which I resample over time to keep the size of the FIFO some target, on average. If my ratio goes too high, my buffer begins to grow and the latency gets large or it overflows, if I am using a circular buffer. If my ratio goes too high, my buffer begins to shrink and eventually underflows. I will periodically find the error between the number of elements in my buffer and the target size. The integral of this error over time, I, will be used to determine R by R=e^-I. The integrator will be set with its initial value I=-ln(R). If the buffer size drifts too high, I goes towards zero and R drops, shrinking the buffer. If it drifts too low, the opposite occurs. What do you think? Should I expect anything weird (oscillation through my feedback path, etc.)?


texts: Formula for calculating IIR filter coefficients for a high order low pass?

texts: I'm trying to make a filter for use in real-time audio processing and I'm trying to figure out how to produce coefficients for a low pass with a steep attenuation curve. I've found a few examples of b0, b1, b2, a1, a2 but I'd like to have the option of a high order filter, which to my knowledge means more coefficients. How are the extras calculated? I'd be extremely grateful for any help on this. Thanks!


texts: Frequency Estimation Using Multiple Signal Classification Algorithm (MUSIC)

texts:
I'm working on side research that deals with signals and data in the time domain. currently, I want to use the multiple signal classification algorithm (MUSIC) which is mainly used for frequency estimation (I've found a lot of resources for DOA MUSIC which is not my interest). The formula of the equation is as follows:-
$$P_x{}(f) = \frac{1}{\sum\limits_{k=p+1}^{N}|s^T{}(f).v_i|^2}$$
where:-
S(f) is a complex sinusoidal vector which is represented as follows
$$e^{j2\pi kf}$$
T is the transpose of the complex vector
Vi as proposed in the book, it represents the eigenvector in the noise subspace
f is the frequency
It's the same equation in Wikipedia and the book below Wikipedia link
My main issue is I don't know what frequency values should be used in the equation (in fact I'm interested in a frequency range of the output but I don't know how it should be evaluated in the equation).
My second issue is how that could be done mathematically especially when I do the dot product between the frequency vector s(f) and the eigenvector Vi since what I've made is I've created a frequency vector of size n (same as eigenvector size) and I made the dot product between them, I'm not sure if I just need to do the dot product for a single frequency vector or all the possible values of k.
To summarize things up, I have a list of values that represents the signal i.e. 300 elements in a given time i.e. 30 seconds. and here's what I've did

I've created an autocorrelation matrix using the Toeplitz matrix algorithm.

I've computed the eigenvalues and eigenvectors of the autocorrelation matrix.

I've sorted the eigenvalues in descending order and then the eigenvectors based on eigenvalues.

I've decided the value of P, as it's described in the book in which there should be a noticed difference (huge) in the eigenvalues.

I've splitted the eigenvectors into principal eigenvector and noise eigenvector (we're only interested in the noise eigenvector).

I've created a frequency vector based on the equation for all the possible values of K for a given frequency.

I've made a dot product between the frequency vector and the noise eigenvector which have the same length.

Here's a pseudo-code of my implementation, but I don't think that it's correct since I don't understand the equation mathematically well:-
N <- the number of samples in the signal
P <- the inedx of the value which splits the singal into principal sinal and noise signal
SRate <- the number of samples per second in our case fps
Vi <- eigen vectors in the noise subspace
s(f) <- the complex sinusoidal vectors (frequency vector)
s(f)^T <- the transposed matrix

sig --> signal

acm <- autocorrelation matrix of sig
eval ,evec <- compute the eigen values and eigen vecotrs of acm

evec <- sort_decsending_based_on_evals(eval, evec)

v <- evec[p:N]

slices <- SRate * N

for f_i <- 0 to slices:
    f = f_i / slices
    ω = 2 * π * f
    // constructing the frequency vector
    for i <- 0 to N:
        // e^(j*omega) was computed using euler's formula
        s[i] <- cos (i*ω) + j sin(i*ω)

    sum <- 0
    for k <- p + 1 to N:
        sum += | s^T . v[k] | ^ 2

    res{f} <- 1/sum

This is the reference of the book in which the algorithm exists page 91 :
book


texts:
I am generating noise in two different cases with the same SNR once using AWGN command which is a built in command and second time using randn function with the correct scaling however I am getting different noise power in each case.

Following is the script to generate noise using different techniques. Can anyone help me why I am not getting the correct result?

Method 1

  A=[1-1i 1+1i 2-1i 2+1i  3+1i  3-3i]
    Norm=sum(abs(A).^2)/length(A);
    flipud(A);
    snrdb=0;
    snrlin=10^((snrdb)/10);
Noise=sqrt(Norm/(2*snrlin))*(randn(size(A,1),size(A,2))+1i*randn(size(A,1),size(A,2)));

Method 2

B=awgn(A,snrdb,'measured');
Noise2=B-A;

The result is as follow:

>> sum(abs(C).^2)/length(C)

ans =

    3.6551

>> sum(abs(Noise).^2)/length(Noise)

ans =

    6.4898

Power of signal is

>> sum(abs(A).^2)/length(A)

ans =7

This seems that randn does better job in estimating the power of signal and awgn this case but in different cases the results are varying and might get different result. Can anyone tell me how I can change my code that both noise generate using built in command and randn gives me the same result.


texts:
I am trying to figure out what seems very trivial thing about down sampling and up sampling. The problem that I am facing is how to find the effects of up sample on a signal followed by a filter and/or a filter followed by down sampling. Here is a problem related to this from the Oppenhiem and Schafer book.

I am trying to see what happens to x[n] after it pass the hlp[n], considering that hlp[n] has wc = pi/4. Will there be a gain on this or just the frequency is affected?

Please help, Thanks in advance


texts: BER vs Eb/N0 for BPSK signal with depth of modulation 120 degrees (+-60 degrees)

texts:
I want to draw graph of function BER vs Eb/N0 for BPSK modulation with depth of modulation is from –60 to 60 degrees (or from 0 to 120 degrees). Can I use for that MatLab's tool 'bertool', or I need to use something else?

Thank you!


texts:
You did it using the norm of the function(signal).

I want to measure this symmetry by subtracting each point of my signal from the flipped one. My question is about the way to normalize this difference. Can standard deviation of the signal be a good option? ( I do not want to divide it by the value of max).


texts: How can I automatically identify and remove a fundamental frequency and all of its harmonics?

texts:
I have some microscopy data that is contaminated by a heartbeat artifact that I'd like to remove. The data consists of a large time series of images captured at ~60Hz.
Here's a small example clip in GIF form:

I have taken the average pixel intensity over time, and computed the periodogram using Welch's method:

As you can see there is a sharp peak at ~1.8Hz which is likely to correspond to the heart rate (~108 beats/min). There are also a bunch of harmonic peaks at integer multiples of 1.8Hz. The exact heart rate is likely to vary from dataset to dataset, but I can specify a biologically plausible range as shown by the shaded area on the periodogram.
What I'd like to be able to do is:

Automatically detect the fundamental frequency corresponding to the heartbeat, and all of its harmonics
Filter the data so as to remove the fundamental and all harmonics.

At the moment I can solve point 1 very crudely by finding the largest peak in the periodogram, then multiplying it by $1, 2, ..., N$ where $N$ is the estimated number of harmonic peaks, but I'm sure that there must be a better method than this hack.
Regarding point 2, I came across this question which mentions using a comb filter to remove a fundamental and all of its harmonics. Is this the best method to use? One important consideration is that I will have to apply the filter to each pixel timeseries in a large array, so a computationally efficient method would be highly desirable.
Example data

I've uploaded a 1D timeseries of example mean
pixel values in .csv format here.
I've also added a sample of 1000 frames (spatially decimated 2:1) in a compressed .mat file here


texts: non-uniform sampling on the unit circle - applications

texts: If there is an N-point sequence, it has a Fourier Transform.  DFT of this sequence is the sampling of the Fourier Transform at N equally spaced points on the unit circle.  What happens if we take some M non-uniformly spaced samples on the unit circle with M < N, M = N and M > N.  Any studies done on this.  If yes, what are the applications?


texts: I have a 10-bit ADC with 8.5 bit ENOB. I have 2.5 quantization levels of thermal noise in front of the ADC. Is there a reason to process more than 8 bits from this ADC?


texts: Low-pass filtering bounces when signal doesn't start at 0

texts:
Say I have an array representing a signal, going from 30 to -30 in 1000 increments, with a little bit of noise on those values. If I apply any king of low pass filtering (Besset, Butterworth, Savitsky-Golay, etc), then the 1st point is now very near zero and the next few points reach progressively for +30 after some wild oscillations depending on the filter and its parameters.

It's been a while (ahem, decades) since I last delved into the maths of DSP, so I'm a bit rusty and looking to know what I should do about it. I thought that the solution was to window my data before filtering it, but it doesn't seem to change much.

Obviously what I want to end up with is a nice 30 to -30 ramp, with the noise removed.


texts: Filter Percentage Uniform Noise from a DC-signal?

texts:
I'm not good with signal processing but i've looked around and have got no clue how to approach this.

My Question is, If there is a Static value present - that is being corrupted by percentage uniform noise, is there a way to estimate the original static value given N corrupted measurements.

i.e

For a Constant X,
On each measurement, K = X + X*N
Where N is a noise signal with a Uniform PDF between bounds +0.50 to -0.50.

How to approach processing this noisy signal to estimate the original DC-value ?


texts: Calculating filter parameters from system response

texts:
In trying to answer another question, I got stuck applying what I thought would be a valid way to answer it: using the system function. Here's what I did:

First, I put the question

Y[i] = Y[i-1] + ALPHA * ( X[i] + Y[i-1] )

into DF I:

Y[i] = ALPHA * X[i] - (ALPHA-1) * Y[i-1]

This tells me that the system function is:

H(z) = ALPHA / ( 1 + (ALPHA-1) z^-1 )

Substituting:

z=e^(jw)

and

|H(z)| = .5 (or .707)

should allow me to calculate ALPHA, but I only get complex values for ALPHA. What am I doing wrong?


texts: Creating a 4th order Butterworth filter R

texts:
I have sampled my data at 2000 Hz and want to generate a fourth order low-pass zero-lag filter with a cut-off frequency of 150 Hz.
I have put the code below that I have created but I have the following questions. Given that I have used the filtfilt function do I need to change the number to “2” so that it is a second-order and should I also double the cut-off frequency i.e change “W” to 0.3 from 0.15.
‘’’rawdata = data$Moment unfiltered
bf = butter(4,0.15,type = c("low"), plane = c("z"))
filtereddata = filtfilt(bf, data$Moment unfiltered)’’’
Thank you in advance!


texts: Convolution for audio, how to control amplitude?

texts:
I am implementation a convolution reverb method for audios array in Python using scipy.signal.convolve. I choose the fft mode, which essentially a product of the rfft of the two signals then ifft back.
What I am not sure is how to perform the convolution without changing the amplitude of the original signal. My audio sample is normalized to -1. to 1.0. And the impulse response (ir) is not normalized but has a peak value near 1.0.
A direct convolve(signal, ir) will result a new signal that has the amplitude of over -6.0 to 6.0. This of course results in a massive clipping during playback. Of course I can reduce the amplitude of the ir to bring down the level. However, I don't know what is the mathematical relationship here in order to achieve the same amplitude.
Thank you in advance for your advice.


texts: OFDM: How to do the carrier recovery of the quadrature(I/Q) demodulation

texts:
For an OFDM system, after the OFDM modulation, the quadrature modulation is used to move the signal from the baseband to the pass-band.
But when it comes to the receivers, I can't find any way to make the carrier recovery, not only the frequency but also the phase.


texts: What is the difference between sampling over $-4\pi$ to $4\pi$ and 0 to $8\pi$

texts:
I see no difference in magnitude, but the phase plots vary significantly while using NumPy and MatPlotLib to plot.
Looking at the phase plots, I feel like I'm missing something important. Also, I'm a relative newcomer to the field of DSP, so if it's something obvious, even then please explain it as an answer rather than a comment.
My doubt is whether there is some error in the phase that is being introduced by NumPy while taking the FFT of a real-valued function or if it is expected to be this way.
Plot using 512 samples over $[-4\pi,4\pi)$:
$[-4\pi,4\pi)$)" />
Plot using 512 samples over $[0,8\pi)$:
$[0,8\pi)$" />
This is my code for generating each of these DFTs of the amplitude modulated cosine wave:
t1 = np.linspace(-4*np.pi,4*np.pi,513)[:-1]
y1 = (1 + 0.1*np.cos(t1))*np.cos(10*t1)
Y1 = np.fft.fftshift(np.fft.fft(y1))/512

t2 = np.linspace(0,8*np.pi,513)[:-1]
y2 = (1 + 0.1*np.cos(t2))*np.cos(10*t2)
Y2 = np.fft.fftshift(np.fft.fft(y2))/512

Do let me know if the rest of my code would also be necessary. It's honestly just plot functions other than these.
Also, please do link me and mark it as duplicate if this is one. I tried searching on Google, but I wasn't able to find anything explaining this concept.


texts: Is an image a random variable or a random process?

texts:
Is an image a random variable, where each pixel is a realization of the same random variable? Or is an image a collection of random variables (if an image is of size mxn, there are mxn distinct random variables)?

In the latter case, are the RVs (pixels) considered i.i.d. (even though they are not)?


texts: Find a stable transfer function $G(z)$ such that $|G(z)| = |H(z)|$

texts:
Consider the following causal IIR transfer function:
  $$ H(z) = \frac{2z^3 - 4z^2 + 9}{(z-3)(z^2+z+0.5)} $$

  Is $H(z)$ a stable function? If it is not stable, find a stable transfer function $G(z)$ such that $|G(z)| = |H(z)|$.

Now obviously the transfer function is not stable due to the pole at $z = 3$.
And I would have just thought that the solution was to flip this pole inside the unit circle, so $z = 1/3$.

And so  $G(z)$ should be $  = \frac{2z^3 - 4z^2 + 9}{(z-1/3)(z^2+z+0.5)}$ right?

The solutions however go $  G(z)= \frac{2z^3 - 4z^2 + 9}{(1-3z)(z^2+z+0.5)}$, which seems fine, except for the next step, where they go $G(z) = H(z)A(z)$ where $A(z) = \frac{1-3z}{z-3}$

That part makes no sense to me, since $$H(z)A(z) =  \frac{(1-3z)(2z^3 - 4z^2 + 9)}{(z-3)^2(z^2+z+0.5)} \neq G(z)$$

  Can someone please explain a) what am I doing wrong; and b) why we introduce $A(z)$?


texts: Could you tell me why wavelets without scaling functions can't be used for two-dimensional analysis? What is the role of scaling function?


texts: what does frequencies in non periodic signals mean?

texts:
What do the frequencies in the a Fourier transform of a non-periodic signal mean physically?
Are there another definition of frequency that doesn't include the FT?


texts: extracting MPEG-7 low-level audio descriptors using essentia

texts:
According to https://www.ntnu.edu/documents/1001201110/1266017954/DAFx-15_submission_43_v2.pdf, the audio processing toolbox Essentia includes the 100% of MPEG-7 low-level audio descriptors. However I have been trying to use it to extract them from a .wav file, and I can't find their names in Essentia's documentation.

My question is: Does Essentia actually include those descriptors out of the box (I mean, as a function, as it does with the MFCCs for instance), but maybe with different names than the originals? Or is it that it only includes "general tools" (as Spectrum extractors) to help with the actual extraction of those features?

In either case, could someone please explain here (or refer to a place where it is explained) how to get the seventeen MPEG-7 low-level descriptors USING ESSENTIA?


texts: How to get rid of discontinuities at the start of a signal caused by the filter

texts: I'm processing an audio signal using an overlap add method.Everything at the block boundaries is fine. However I have a discontinuity at the very start of the first block and I'm not sure how to (or if I can) get rid of this.Zero padding the start seems to make no difference as the discontinuity is obviously "triggered" when the filter first encounters the non-zero signal.


texts:
Suppose we have a vector of length N of DFT, which every odd index (in frequency) equals 0:

$$X[k] = 0\quad \forall k:\, k\mod2\equiv1$$

What does that mean for the series in the time domain?


texts: Power Spectral Density vs Frequency Domain separation

texts:
I am just starting to learn about Power Spectral Density (PSD). What I'm reading is that it is possible to separate PSD masks if they are isolated, even if the frequencies do overlap. My understanding is that this is how telephone lines for example can carry both telephone and adsl signals simultaneously.
Now, I also came across a claim that frequency isolation is a "sufficient but not necessary condition" for signals to be separatable. Does this mean that if the frequencies are isolated, then the PSD signals do not need to be isolated? In other words, isolation of one or the other is sufficient?
Please excuse any wrong terminology, and I welcome corrections. I'm quite new is this field, so appreciate any help.


texts: How best to analyse a 130volt +/- 20volt signal?

texts:
I want to measure and log a relatively high variable voltage signal, with minimum impact on the original signal.

Normally the signal should average 130volts and vary by +/- 20volts at 10Hz.  However under different circumstances all these parameters are likely to change.

I would like to log these changes, but the problem is my data logger can only support volatages up to 60v DC or 42v AC.

I have considered using a small transformer, but fear that the relatively slow frequency may not induce a representative current in the secondary coil.  And/or the variable frequency may result in back emf, corrupting the original signal.


texts: How do I decode(demodulation) and get original data from audio sample data

texts:
Here I'm talking about the less-than-ideal situations, in other words, the audio
sample data may have burrs or error bits in it.

Encode and decode theory:
I am using FSK to encode my original binary data (0,1 and lead-byte), and transform the square wave it generated to sine wave. Then I generate a stream of audio samples through mic (for example, a stream of 44100 16-bit integers for every second etc.).
When I decode (c language), the key factor is, 0:1:lead-byte = 1:2:4(or other ratio).

For example:

original binary data: 0|1...0|1.

encoded wave data: lead-byte...lead-byte 0|1...0|1 lead-byte.

The typical ideal samples and bad samples waveform are as below.

good wave

bad wave 1

Now my decode method is simple (just count the samples up and below zero, and get stable lead-byte width first, then use 1:2:4 to calculate 0 and 1's width). Should I use more complex mathematical analysis in decoding? If so, what theory can I use, Fourier? Gaussian Noise? Rayleigh Fading? Goertzel algorithm?

I have to point out that I am totally new in audio stuff.
Any help would be much appreciated!


texts: Random Process at a particular time instance

texts:
I was studying Random Process and I thought I understood what it was all about until I came across this example. Consider a random experiment of tossing a coin with sample space S = {H, T} The sample functions are given by:

$X(t, H) = x_1(t) = \sin(\omega_1t) \tag{1}$

$X(t, T) = x_2(t) = \sin(\omega_2t) \tag{2}$

Now, I understand that random process at any particular time instance say $t_k$ gives me a random variable. Also, the sum of probabilities of a random variable should be 1 i.e.
$$X(t_k,H) + X(t_k,T) = 1 \tag{3}$$
subsequently in general no matter where I sample the random process I should have,
$$\sin(\omega_1t) + \sin(\omega_2t) = 1 \tag{4}$$

But the Equation 4 won't be true for all values of t.

Question 1. Where am I going wrong above?

Question 2. If I take the value of $X(t,\lambda)$ at a particular time $t_k$and at a particular $\lambda_i$ the value $X(t_k,\lambda_i)$ will be a number. This number is equal to the probability of $X(t,\lambda)$ taking value $\lambda_i$ at time instance $t_k$. Am I correct?


texts: What's the Intuitive Description of Circular Symmetric Complex Zero Mean White Gaussian Noise?

texts:
Suppose we have an added discrete noise signal defined as:

$$y[n]=x[n]+w[n], $$ where $w$ is zero mean white Gaussian noise.

Question:

When we say white noise, is it sufficient to say that it is always be uncorrelated with signal $x$, no matter what the distribution of signal is? (I think yes)
Can we conclude that both $x$ and $w$ are independent IF noise is white Gaussian  BUT the signal have unknown distribution ? (I am Confused about it)
Lastly, What the information we glean-out of the model (intuitively) if we are informed about noise that besides being zero mean, white and Gaussian, it is circular symmetric complex as well? (I do know white noise has impulse auto-correlation but confused with circular-symmetric term).


texts: Difference between band pass vs. low pass + high pass?

texts:
I am reading a paper in which they filter EEG signal using a notch filter to remove the power supply component of 50 Hz, then a high pass filter at 0.3 Hz to remove artifacts resulting from breathing, and finally a low pass filter with a cut-off frequency of 40 Hz.

Is there any difference to simply using a band pass filter of 0.3 – 40 Hz?

Additionally, why the notch at 50 Hz when everything 40+ Hz is thrown away, anyway?


texts: How exactly does sound get encoded by a cochlear implant processor?

texts:
I've been wanting to know more about how sound is processed by a cochlear implant processor, how they are programmed, and how they stimulate the magnets and electrodes inside.

How do you take pure tones of varying frequencies, noise signals, and complex speech sounds and turn them into a set of organised electrical signals?

Also, what is meant when people say that cochlear implants can have one channel, two channels, etc? I'm only familiar with mono and stereo channel, but recently I heard that channels refer to the splitting of a single signal into two separate paths. If the signals sound exactly the same in both channels, the result is a mono perception. If there is a difference in signal, the result would be a stereo image.

With cochlear implants, they can have up to eight channels or more. How do these work, and just out of curiosity, what sampling rate and bit depth, as well as bitrate is the sound presented? Is it in Pulse code modulation?

I guess the thing to think about is whether the processor, which is digital converts the analogue sound signals of the microphone into a set of signals that somehow need to be reconverted to another form of Analogue that will stimulate an electromagnet that will stimulate a certain number of electrodes for the brain to perceive sound.


texts:
I have a question about determining time-invariance of a linear system. We are given this system and we need to determine if it is time-invariant or not:

$$y(t)=\int_{-t}^{\infty}x(-3\tau)d\tau$$

Now after doing the shifting $t-t_o$ you get two different outputs. Now that makes sense, but what confuses me is the variable substition that happens in order to check. I'll post a link to the picture on imgur,sorry for this but I don't know how to do it otherwise.

https://imgur.com/a/CBbe4Sp

Basically my question is why is the substittution $-3τ -t0$ and not just $τ -t0$. Also why are we substituing both the Z1 as well as Z2, isn't the goal to make try to make them equal? Like try to modify one to get the same result as the other?.Maybe someone can give me some insight.Sorry for the bad post I need instructions on how to post. Cheers!


texts: How to filter power line noise from ECG?

texts:
currently i'am trying to get rid off power line noise (50/60 Hz) from an ecg.

Specification:

Bandwidth of the ecg 0.05 - 150 (achieving this with a butterworth second order bandpass)
sampling rate 500 Hz

What filters did i try to remove 50/60 Hz noise?

Butterworth bandstop 6th order bandwidth (+/-0.5 Hz up to +/2 Hz)
Notch filter placeing zeros at z = e^(+/-j*2*pi*f_notch/f_sample)
Notch filter with additional poles with same angular as zeros at distance to unit circle with r= 0.5 up to 0.7

What did happen?

Butterworth filter causes ringing in the ecg signal
Notch filter with only zeros causes ringing in the ecg signal
Notch filter with additional poles causes an "overshoot" in the QRS-complex

What is the problem? (I think)

Actually my ecg signal contains information arround 50 Hz due to the phase response (and group delay) of the applied filters, i get this kind of distoration. I mean you clearly remove information from the signal it makes absolutely sense it gets "destroyed".

My question (finally)

Is there are a "better" way to remove 50/60 Hz noise? Secondly i was thinking to take the ringing effect and try to smooth it out (maybe savitzky golay filter?)

Best regards


texts: Goertzel algorithm: Relationship of magnitude

texts:
I've written a quick test app that uses the Goertzel algorithm to determine if a given frequency is present in the signal. This is to pick up DTMF tones and various other signals. The app appears to be working. I'm using the 'magnitude squared' approach in my final stage.

http://www.embedded.com/design/configurable-systems/4024443/The-Goertzel-Algorithm

Now, how does this magnitude relate to all the other variables (especially the amplitude of the source signal)? Most of the implementations determine the presence of a given frequency if the magnitude exceeds some threshold. I want to know how this threshold is determined, and how would I go about selecting a suitable threshold.

EDIT:

OK, I think I've found the answer (I hope). Most implementations of the algorithm state that is will return the 'relative power'. I couldn't find any that stated what it was relative to.

This morning, it dawned on me (pun intended) that it's probably relative to the input signal (duh!). So, if we calculated what the power of the input signal is over all the samples taken, and then compare this to the relative power given by the algorithm, we should have a good idea of the contribution of our target frequency?

Is my understanding correct?


texts: Suggest me some good references for sonification

texts:
I need good reference books/tutorials/doccumentations on sonification.

Please suggest me some books or other references from which I can learn about sonification and algorithms of sonification.

Thanks.


texts: I'm working on a project to stream video between 2fpgas by using QAM. While mapping symbols with different amplitude pulses(PAM) a question popped up in my mind. In real life how does communicating wireless devices know which amplitude levels to agree on for digitally modulating-demodulating bitstreams? Is there a standard for that?


texts: How to Avoid Perspective Distortion Buildup in a Planar Mosaic

texts:
Registering images for map or other planar mosaic from a handheld camera or a drone requires projective transform (homography) as the camera is allowed to rotate between frames. Unfortunately, planar homographies cause mosaic to grow too big with each new image being registered:

Please note I am not trying to make a rotational mosaic, rather a planar mosaic (camera translating) with additional rotation (due to wind, hand rotation etc.).

To avoid this issue, I have allowed homography only for image-to-image registration, but the images are then globally aligned using simple translation, that minimizes distance between feature correspondences (originally matched using homography):

As you can see, this produces big misalignments.

I found that Microsoft ICE (Image Composite Editor) offers "Planar Motion with Perspective" option that tackles this problem:

Here is the result from ICE:

I can't figure out how they solved the problem, i.e. how to combine planar mosaic with perspective distortion and avoid the distortion buildup.

I tried splitting homography transform to in-place perspective distortion without translation, then translating the distorted point:

$$\textbf{x'}=\begin{pmatrix}1 & 0 & t_x \\ 0 & 1 & t_y \\ 0 & 0 & 1\end{pmatrix}\begin{pmatrix}h_{11} & h_{12} & 0 \\ h_{21} & h_{22} & 0 \\ h_{31} & h_{32} & 1\end{pmatrix}\textbf{x}$$

Unfortunately, this still produces the distortion buildup.

Another option is to introduce some error term in the bundle adjustment (optimization) routine that discourages too high values of perspective distortions. The problem with extra error test is that it complicates the model and it is hard to find proper scaling for the additional error factor to properly weight other variables.

Yet another solution may be to keep image positions fixed and adjust only perspective distortion. This however is not a globally optimal solution.

Any ideas?


texts:
I'm having trouble finding information on the topic I'm interested in, and maybe my problem is just how I'm wording things.

So, let's say I have a repeater-based communication system, where one node in the system receives a message, and then passes it along to the next node.

For this example, let's say we have nodes 1, 2, 3, 4, 5, and 6. Also let's assume nodes can always receive the message from their neighbors, but may or may not receive it from further nodes. (ex. 3 will always receive from 2 and 4, and it might receive it from 6)

I'm wanting to find some information on different wants to determine the path.

For instance, i can always assume the path is fixed, where the message goes from Node 1 to 2, 2 to 3, etc. so the path would be 1-2-3-4-5-6. By doing this, the message only needs to contain the node number of the node sending the message

But what if I want to send the path that I want the message to go through, with the smallest amount of data possible? For instance I might want one message to go through the path 1-2-5-6, and another message to go through the path 1-5-2-3-6.

My thoughts to do this is to use a lookup table that is saved on all nodes, but the problem arises when I have a lot of nodes. If I have 100+ nodes, my table is huge and therefore the data I send to associate the path to the table gets large.

Could anyone point me  to the right direction to even look at this subject matter?


texts:
My system is :

A step by step wheel returns its position each 150ms. I want to pilot an horizontal infinite scroll bar with the wheel. The scroll bar needs velocity as input to move.
Moreover the faster the wheel rotates the farthest the scroll must move. That means that if I turn the wheel one turn quickly, the bar must go further if I had turned the wheel one turn slowly.

How can I compute the right velocity to give as input to the scroll bar with sampled position of the wheel ?


texts: what exactly does scale mean in scale-space theory?

texts:
In image processing, scale-space is a technique to represent an image at different scales. But this really confuses me, since I always think that the word scale just means the size of an object, so the scale-space representation of an image shall have different sizes of the same image.

However, it turns out that the scale-space representation is composed of Gaussian-convolved images with different $\sigma$, the image size doesn't change at all, why? And further, how $\sigma$ in the Gaussian function determines the scale of the convolved image ?


texts: What Are Linear and Circular Convolution?

texts:
I have some basic understanding of signals and convolution. As far as I know it shows the similarities of two signals. Could I get some explanation in plain English of:

what are the linear and circular convolution
why they are important
practical situation where they are used


texts: Maximum NCO frequency and filter design for down-conversion of complex input signal

texts: I am trying to perform digital down conversion by using the DDC from Matlab dsp toolbox. My question is that the bandwidth property is limited to Fs/2 and the filters created are real filters. If the input is comlex then the frequency range will be from 0 to Fs and the decimating filters need to be complex but the toolbox does not allow that. Is this some issue with the Matlab toolbox or is this a theoretical limit that DDC NCO frequency cannot exceed Fs/2 even if the input signal is complex?


texts: Bluetooth locating with RSSI and only one receiver

texts:
Is there a method to find the relative direction of a Bluetooth device to the receiver?

A quick Google search on "Bluetooth positioning" returns several methods using triangulation, but I don't want to use multiple receivers. Also, I don't need to know the distance, only the relative angle to the receiver is enough.


texts: Retaining filter state for processing data in chunks

texts:
I am have data in 1 minute chunks but I would like perform filtering on these as if the data was one contiguous signal. How might I hold the filters state between processing data chunks to achieve this.

I am using an FIR in MATLAB with the filter() function.


texts: How do I window my pwelch function on matlab. I only get an average of the first window

texts:
I am implementing a signal processing algorithm to determine tremor and non tremor activity in a signal. The tremors of interest are in the range of 4-8Hz and I need to find the max value for each window but i am only getting one plot and when i use a for loop with a window the same pwelch signal outputs like 6 times. Here is my code:

wind = 2;
ws = wind * samp_freq;
trem_feat = (filt_trem);

for  i = 1:ws:L-ws
    window = i : i+ws;
    [pxx, f] = pwelch(trem_feat ,500 ,300 ,500 ,samp_freq);
    figure;
    plot(f, 10*log10(pxx));
    title('PSD');
    xlabel('Frequency (Hz)'); ylabel('PSD (dB/Hz)');
end

How do i get the average pwelch of all windows?


texts: So if I have an image of the size: 1920×1080 and I want to apply a convolution filter of 3x3. Take the formula for DFT $N^2$ and FFT $N \log N$. What should I fill in for N in this case?


texts: Extracting a DC value from a sum of sinusoids

texts:
I guess it might be a basic questions but here it goes anyways.

How could I extract a DC value from a sum of sinusoids, i.e.:

$$v(t) =\widetilde{v_{dc}} + \sum_{n=1}^{\infty}\sin(\omega_nt+\phi_n)$$

I tried using low pass filters in Matlab or just mean of the signal but it's not quite what I'm looking for. I'd also like to have a mathematical way of expressing this.


texts: Range and motion in received signals by radar

texts:
I've got two signals which are consecutive received signals by a radar system. The signals' values are complex (in form $x+\textrm iy$).
I have generated a trasmitted signal using the info available.

Now,  I have to calculate the range and motion(if any) of the objects using comparision of these signals in MATLAB. Comparing the real part using xcorr in MATLAB, I got the range of the objects.
I am strugling with motion of the objects. How can I find if the objects found (4 in my case) were moving or not. I think it has to do something with the phase comparision of the signals.

Any help is appreciated.


texts: BIBO stable LTI system frequency response for this input signal?

texts:
For an LTI system with bounded input and bounded output, I have the input $$x(t) = 5 + \cos(12t+\pi/4)$$ and output $$y(t) = 6\sin(12t)$$

It is said that the magnitude of the frequency response $|H(j12)| = 6$. I don't understand how this is true. Obviously the second term of $x(t)$ is scaled by $6$, but the first term needs to be scaled by $0$, right? Which means this LTI system isn't possible. Can anyone clarify?


texts:
The Book "Digital Modulation Techniques" by Fuqin Xiong groups digital modulation schemes into two classes: "constant envelope" and "nonconstant envelope."  However, I was surprised by the inclusion of PSK in the "constant envelope" category.  This includes BPSK and QPSK, which are equivalent to PAM and 4-QAM.  This doesn't seem right to me, since both PAM and 4-QAM will have  transitions through zero, unless infinite bandwidth is assumed.

My question then, is what is the definition of a constant envelope modulation?  Is it simply having all symbol constellation points at a fixed energy?  If so, it seems like constant envelope modulation may result in a waveform that is not constant envelope.


texts: MFCC, uncertain of my results and algorithm

texts:
I've started to implement an MFCC algorithm on an embedded DSP, using the CMSIS DSP library. However the results I'm getting don't appear consistent. This is my algorithm so far:

Gather 32ms PCM sample from microphone, 512 samples long. 16000hz
Apply a Hamming window.
Compute the 512 FFT using this.
Compute the magnitudes to get power spectrum.
Generate a set of 30 filter banks - using 32 points on the mel spectrum
Convert my power spectrum magnitudes from 16bit ints to floats.
Get the dot product for each filter with the power spectrum.
Log each dot product.
This is the uncertain part. The DSP library I'm using only has DCT IV, and this only supports sizes 128 and up. My set of logs is 32. So I've opted to use a Real FFT function, the results I get from this vary quite a lot. If I speak the same word in the same tone.

Could this be what is hampering my result? Real FFT provides an output of both real, complex and conjugate data (described in link above) - which would I use for the MFCC? - Should I even expect my MFCCs to be similar or is it expected to vary a lot? My results tend to vary between -20.0 and 20.0.

Thanks in advance.


texts: How to downsample a fourier transformed signal?

texts:
I have a signal of length 100000 timestamps sampled at a frequency of 25kHz. First I apply a high pass filtering at (300Hz) and then do the Fast Fourier Transformation. Then the absolute values are taken from the returned complex values. Finally, the first half of the returned signal (of length 50000), is downsampled to 1000 samples.
signal_chunk_filtered_new = butter_highpass_filter(signal_chunk, 300, 25000, 3)

downsampled = signal.resample([abs(i) for i in fft(signal_chunk_filtered_new)][:int(chunk_size/2)], 1000)

My current implementation uses the scipy resampling implementation for downsampling which uses the Fourier method. https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html
I was simply considering the Fourier transformed signal as another signal which needs to be downsampled. But when I think about it in-depth, it is like I'm applying Fourier transform on top of Fourier transform. Is this method wrong? If so which downsampling method should I use?


texts: Solving Sparse Model with given Dictionary Using LASSO

texts:
I was trying to solve a problem where the basis matrix contains the components of $\sin(nx)$, $\cos(nx)$, $\sinh(nx)$ and $\cosh(nx)$.

Say the $n$ varies from 1 to 100. While solving the lasso linear regression, I should get the sparse coefficients on say $n=10$ for each of the components. To be clear say the data obtained from the equation is $y= \alpha \sin(10x) + \beta \cos(10x) + \gamma \sinh(10x) + \delta \cosh(10x)$.

I am getting some arbitrary coefficients and that is not representing the equation.

What am I doing wrong?


texts: Find response of discrete time LTI system given input and impulse response

texts:
For this question the guide says to use $Y(e^{j\omega}) = H(e^{j\omega}) X(e^{j\omega})$.
I have been able to find the discrete time Fourier transform of the impulse function as $$H(e^{j\omega}) = \frac{1}{1-\frac12 e^{-j\omega}}$$
Finding the same for the impulse function is what I do not understand.


texts: How does the Hilbert Transform extract the in-phase and quadrature phase components of a real passband signal?

texts: I am taking an intro to radar signal processing course and we are discussing the Hilbert Transform and how it take a real input and its output is analytically complex. I'm trying ) how a radar recextracts $p(t)$ and can $q(t)$ frathe in-phase and quadrature phase components ssband signal $t$, (sayo underst)and how and can give insight into amplitude and phase. Could someone shed light on this topic, or perhaps provide a block diagram of how this process works?


texts: Expected Value of a sequence with two random variables

texts:
If I have a signal of the form $x\left(n\right)=Acos\left(nω+ϕ\right)$
where
$\omega \in \left[\omega -\lambda ,\omega \:+\lambda \:\right]$ is a uniform random variable
and
$\phi $ is also a uniform random variable $\left[-\pi ,\pi \right]$.
I need to find the mean and the autocorrelation.
So we know the mean $E\left[x\left(n\right)\right]\:=\:AE\left[cos\left(n\omega +\phi \right)\right]$.
I think the $E\left[x\left(n\right)\right]\:=0$ regardless of the random variables just because of the periodicity of the cosine sequence.
Now I would also like to find the auto-correlation function
$R\left(n,k\right)=\:E\left[x\left(n\right)x\left(n+k\right)\right]=A^2E\left[cos\left(n\omega +\phi \right)cos\left(\left(n+k\right)\omega \:+\phi \:\right)\right]=\frac{A}{2}^2\left\{E\left[cos\left(2\left(n\omega +\phi \right)+k\omega \right)\right]+E\left[cos\left(k\omega \right)\right]\right\}=\frac{A}{2}^2cos\left(k\omega \right)$
I am not sure if am solving this the right way, can somebody please tell me if am doing something wrong? because I think the random variables must have a role in this


texts:
To understand the idea of sparse signal recovery described in Sparsity and Incoherence in Compressive Sampling, I decided to build a toy problem. Suppose the sparse signal that we want to reconstruct $x^0$ is $n=1000$ dimensional and only $S=9$ of its entries are nonzero. For simplicity, I pick the transform matrix $U$ to be a DFT matrix and I randomly sample $m$ values from $Ux^0$ (transform coefficients) to use them as our measurements $y$. In other words,
$$y = U_{\Omega} x^0.$$
Based on the theorem given in the paper, signal recovery is possible with constrained $\ell_1$ norm minimization if the number of measurements $m$ satisfies the following condition:
$$m \geq C \mu^2(U) S \log n$$ for a positive constant $C$. In my toy problem, I know that $\mu(U)=1$, $S=9$ and $n=1000$, but I do not know the constant $C$! I looked at the paper over and over again and I think I am missing something. Could you please explain to me how to determine the value of $C$?


texts:
I have a stream of music that I'm reading(from a microphone, from an online streaming source or whatever). And I want to have a system that will detect when a song changes.

The question is, is this possible with any kind of accuracy?(even something 75% rate of detection would be helpful).

Alternatively, would it be possible to detect large enough changes in the music(so it might be a song change but it could just be a switch from a verse to chorus or chorus to bridge, etc...).

I am not asking for an exact algorithm -- though obviously if you know one then I obviously wouldn't mind -- but even a general approach on how to tackle such a problem would be helpful.

All pointers welcome, scientific papers, open source projects, some random guy's blog, whatever.


texts: Infinite extent of spectrum, but also in time in Oppenheim's Discrete Time Signal Processing?

texts: In Oppenheim's Discrete Time Signal Processing there's on p. 323 no limited band in both time and frequency - wouldn't that violate the Heisenberg Principle?


texts: Complex Filter design using Matlab's FDAtool

texts:
What functions does matlab's FDAtool use to design filters from pole zero placement?

Functions such as zp2sos require the the poles and zeros are either real or complex conjugates, FDAtool does not have these constraints and so can design complex filters.

Is there then a function that given zeros and poles that are not complex conjugates will return the coefficients for a sos iir filter?


texts: Minimum Mean Square Estimator - Equivalent Expressions to Minimize

texts:
Given $ M \in \mathbb{R}^{N \times N} $ which is a Positive Definite Matrix.
Let $ \hat{x} $ the MMSE of $ x $ given $ z $, namely $ \hat{x} = \mathbb{E} \left[ x \mid z \right] $.

Prove the equivenalce of the following expressions:

$ \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} \left( \hat{x} - x \right) \mid z \right] $.
$ \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} M \left( \hat{x} - x \right) \mid z \right] $.
$ \arg \min_{\hat{x}} \operatorname{Tr} \left( M \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} \left( \hat{x} - x \right) \mid z \right] \right) $

Equivalence means all are actually minimized by $ \hat{x} = \mathbb{E} \left[ x \mid z \right] $.

The equivalence of 2 and 3 is easy using the Cyclic Property of the $ \operatorname{Tr} \left( \cdot \right) $ operator.
Yet showing the invertible transformation keeps the solution (Equivalence of 1 and 2) isn't trivial.


texts: Filtering - multiplication in frequency domain

texts:
I am trying to create a simple low pass filter, but I have got what I consider a surprising result when looking at the frequency response of a simple Butterworth filter.

I have copied much of the example below from this other post. I have added some code at the bottom of the script to compare the input and output spectra against the frequency response of the filter. I would expect that the output spectrum $\mathbf B$ should be the product of the input spectrum $\mathbf A$ and the frequency response $\mathbf H$:
$$
\mathbf B = \mathbf H \mathbf A
$$

However, the plot on below shows that the filter actually increases some low frequency components - see how the red line is above the green below around $4\textrm{ Hz}$.

Could anyone please explain why this is?

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
from scipy.fftpack import fft as fft
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

# Filter requirements.
order = 6
fs = 30.0       # sample rate, Hz
cutoff = 3.667  # desired cutoff frequency of the filter, Hz

# Get the filter coefficients so we can check its frequency response.
b, a = butter_lowpass(cutoff, fs, order)

# Plot the frequency response.
w, h = freqz(b, a, worN=8000)
plt.subplot(2, 1, 1)
plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
plt.plot(cutoff, 0.5*np.sqrt(2), 'ko')
plt.axvline(cutoff, color='k')
plt.xlim(0, 0.5*fs)
plt.title("Lowpass Filter Frequency Response")
plt.xlabel('Frequency [Hz]')
plt.grid()

# Demonstrate the use of the filter.
# First make some data to be filtered.
T = 5.0         # seconds
n = int(T * fs) # total number of samples
t = np.linspace(0, T, n, endpoint=False)
# "Noisy" data.  We want to recover the 1.2 Hz signal from this.
data = np.sin(1.2*2*np.pi*t) + 1.5*np.cos(9*2*np.pi*t) + 0.5*np.sin(12.0*2*np.pi*t)

# Filter the data, and plot both the original and filtered signals.
y = butter_lowpass_filter(data, cutoff, fs, order)

plt.subplot(2, 1, 2)
plt.plot(t, data, 'b-', label='data')
plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend()

plt.subplots_adjust(hspace=0.35)
plt.show()

def calculateFFT(time,signal):
    N=len(signal)
    df=1/((time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [2.0/N*abs(i) for i in fft(signal,N)[0:N/2.0] ]
    return frequencies,fftValues

plt.subplot(2, 1, 1)
originalfreqs,originalFFT=calculateFFT(t,data)
plt.plot(originalfreqs,originalFFT,"g",label="original")

filteredfreqs,filteredFFT=calculateFFT(t,y)
plt.plot(filteredfreqs,filteredFFT,"r",label="filtered")
plt.legend()


texts: Moving average filter using DF-I and DF-II

texts:
I understand that the moving average filter is simply the average of a number of points from the input signal to produce the output signal. If x[] is the input and y[] is the output and M is the number of points, then the filter is just $$y[i]=\frac{1}{M}\sum_{j=0}^{M-1}x\left[ i+j \right]$$
This can be done through many ways, but I want to ask about using Direct Form I and Direct Form II. How do I follow the DF-I and DF-II structures to create a moving average filter? Don't I need a two separate multipliers like $$b_{M}$$ and $$a_{N}$$ to use the DF-I and DF-II, which moving average filters don't necessarily have? How about if I want to use DF-I and DF-II to create an autoregressive filter? Isn't it even more incompatible with DF-I and DF-II?


texts: Help in problem formulation for estimation of image as a feature vector - SISO or MIMO FIR channel model?

texts:
Based on the paper Blind Image deconvolution:

A feature vector is a list of numbers used to represent an image. The feature vector for my case takes values as symbols $-1,1$.

An instance or an example of an image is represented as a feature vector of $d$ feature values $\mathbf{x} = [x_1,x_2,...,x_d]$. I have $N$ examples
$$\{\mathbf{x}_i\}_{i=1}^N$$ in the form of a database. So, the database contains $N$ rows of feature vectors each of $d$ length.  One such pattern (feature vector) is
$$\mathbf{x} =  [1, -1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, -1, -1, -1, 1]'$$
where $d = 30$
The database thus has $N$ rows of examples each of $d$ length :
$$ \mathbf{x_1} = x_{11},x_{12},....,x_{1d}$$
$$ \mathbf{x_2} = x_{21},x_{22},....,x_{2d}$$
$$ ::$$
$$ \mathbf{x_N} = x_{N1},x_{N2},....,x_{Nd}$$

I don't know how to use the above as input to the channel. In general,
the output is,
$$\mathbf{y} = H^T x$$

Assuming that the impulse response is that of FIR, then I am facing difficulty in how to represent the channel. Should the source input be each component of the feature vector or an entire example?

if the input in the form of the feature vector is transmitted via a channel whose impulse response is modeled as moving average (finite response, FIR) then inputs and how many outputs should there be? The output would be the convolution of the input and the impulse response. Would the channel be Single input-Single Output FIR or Multiple Input Multiple Output FIR ? The objective is to estimate the channel parameters and the input using estimation methods from the output observation only. Would the estimation be performed using one example or all the examples?
Should there be $d$ channel coefficients?
A mathematical representation of the system will really help to clear the concept. Thank you


texts: Denoising Oscilloscope data the right way [python3]

texts:
I need to plot large oscilloscope data sets. They are too large (100 million Values) to plot all of them in a for loop. If I try to do that it results in a MemoryError. So with the help of Stackoverflow I came up with a solution that isn't pretty but it works.

Original Dataset:

To reduce the amount of data to be plotted I can remove all rows of the dataset, that have a y-value between 1 and -1. This can be done comfortably with pandas:

df is the dataset in form of a pandas.DataFrame().

df[1] is the column which is to be examined.

df.drop(df[(df[1] < 1) & (df[1] > -1)].index, inplace = True)

This results in:

This reduces the length of the dataset from 100 million to about 200 thousand. My Question is: Is there a way to reduce data length and keep its original look (prevent the ugly cut-off)?

This stack seems to work differently than Stackoverflow. Pardon me for all wrongdoings.


texts: Electret microphone RC filter design contradiction among information sources

texts:
I'm building a circuit for an electret microphone and I want to build a bandpass filter around the op-amp.
I was using this source until I have found that most sources ( e.g. here and here) indicate that RC filter needs to be grounded.
I can't think of a way check this (FFT doesn't seem to indicate whether the first source is correct or not). Any ideas on how to check or information which source is correct?


texts: Sufficient conditions for exact signal recovery using OMP?

texts:
For a compressive sensing model :
$$y_{_{MXN}}=A_{_{MXN}}x_{_{NX1}}$$
where $x$ is $K$ sparse, what is the sufficient condition for Orthogonal matching Pursuit (OMP) to exactly recover the data for both noisy and noiseless case ?

There are multiple statements, like

$M>Klog(N-K)$ ;
$M>Klog(N/K)$ ;
$M>2Klog(N-K)$ ;
somewhere it's $M>4Klog(N-K)$.

Should there be some constant ratio between $M$ and $N$?


texts: What do these values in the freeverb diagram mean?

texts:
I am implementing a reverb based on the Freeverb implementation detailed here. There is also a more in-depth description of the LBCF in question here, but I am not able to find the information I want at that link.

In the diagram given, they use a lowpass-feedback comb-filter, denoted by LBCF(f,d,N). N is clearly the number of samples. I suspect d means the feedback gain, though I'm not sure if a value of .2 makes sense for that. I have no idea what f means.

So here are my questions:

What does d mean?
What does f mean?
What cutoff frequency is used for the lowpass filters?


texts:
I have some images which have some objects of my interest. The objects have very low contrast w.r.t. background and I would like to detect the edges of these objects in these low-contrast images. Any suggestions?

An example image is shown below (objects are the rectangles in the right side of the image and the edges are their boundaries).

Thanks.


texts: Time-frequency analysis of non-sinusoidal periodic signals

texts: Given the history of the sum of a time-varying mixture of periodic signals, say square waves, how would you efficiently estimate the number and frequencies of components active at a particular time? The amplitudes and frequencies of the components are arbitrary but fixed real numbers; if a component is active at a certain moment, it will retain its amplitude and frequency if it is active again.


texts: Why will the power of signal become larger in the receiver ,but not smaller? the channel is a rayleigh fading channel

texts:
$\mathbf H_{AB}$ is a rayleigh fading ,$7$ is the antenna number for transmitting in the transmitter A, and $2$ is the antenna number for receiving in the receiver B

  H_{AB}=sqrt(1/2)*(randn(2,7), j * randn(2,7))

$\vec f^u_A$ is a $7$ by $1$ precoding direction,$(\vec f_A^{\ u})^H \vec f^u_A=1$

$p_A$ is a transmitted power of $\vec f_A$,assume $p_A=2.5W$,that is $\sqrt{p_A} \vec f^u_A=\vec f_A$

Now i sent the message from A to B,that is the SNR of receiver is

$SNR=\frac{| H_{AB} \  \vec f_A|^2}{N_0}=\frac{\vec f_A^ {\ H} H^H_{AB} \  \vec f_A \ H_{AB}} {N_0}=\frac{10W}{1W}=10$

I want to ask why will the received power be larger than the transmitted power ?

I think the power should decrease after going through the fading channel.where do i wrong?


texts: Frequency change detection in step-wise frequency chirp

texts:
I am driving a system with a step-wise frequency chirp from 50-400 hz with a frequency step of 0.01 hz. I am having difficulty identifying the points of frequency change as the data collected has an indeterminate amount of time at 50 Hz and an indeterminate amount of time at 400 Hz. Since the time between frequency changes is 0.1 seconds an fft calculation doesn't provide the proper frequency resolution to determine this change. I also tried calculating the gradient of the frequency chirp with inconclusive results as well as attempting a lock-in calculation. I wrote some sample code that has a randomly located frequency step followed by another 0.1 secs later that I have been trying to diagnose with varying success. Ultimately I am trying to identify the beginning and end of the frequency chirp since at that point the step changes are known but this change point identification was attempted to identify the point of change between 50 and 50.01, then backtracking to find the beginning of the series. Any suggestions would be very much appreciated.

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

fs = 10000
t = np.linspace(0,1,1*fs)
x = np.zeros(len(t))

step_loc = int(2500 + 5000*np.random.random(1))

x[:step_loc] = np.sin(50*2*pi*t[:step_loc])
x[step_loc:step_loc+1000] = np.sin(50.01*2*pi*t[step_loc:step_loc+1000])
x[step_loc+1000:] = np.sin(50.02*2*pi*t[step_loc+1000:])

grad = np.gradient(x,t[1])

lock = x*(np.cos(50*2*pi*t) + 1j*np.sin(50*2*pi*t))
b, a = signal.butter(5,30/10000, 'low')
lock_filt = signal.filtfilt(b, a, lock)

plt.plot(t,x)
plt.plot(t,grad*3e-3)     #this is offering very inconsistent identifiable results
plt.plot(t,lock_filt*1e1) #although this can identify the general area of the random frequency change i
                      # dont know how to extend it to identifying the actual point


texts: AutoCorrelation using FFT of power spectrum

texts:
I have a chunk of data (samples) that need to do autocorrelation using FFT method. All data are real and lies between -1 and +1. What I have done follows:

Zero pad the window with the length equals the samples' length.
FFT of result from 1
Result from 2 has Real and Imaginary parts. I used them to find Power Spectral Density by taking Re*Re + Im*Im.
Take the inverse Fast Fourier Transform.

Now here is my question: The result from part 4 contains Real and Imaginary parts. Is the magnitude of it equals the autocorrelation result? If so, since my input contains positive and negative values, the magnitude of Result from part 4 means the autocorrelation result cannot be negative, which doesn't make sense.

Or should I take only the Real part of result from Part 4 as the autocorrelation result? My input is real and so the autocorrelation shouldn't have any imaginary part, right?


texts: MATLAB IIRs -- extracting and scaling the coefficients of second order sections

texts:
suppose I generate the following filter:

Fs = 32000;  % Sampling Frequency
N     = 4;    % Order
Fpass = 400;  % Passband Frequency
Apass = 1;    % Passband Ripple (dB)
h  = fdesign.lowpass('N,Fp,Ap', N, Fpass, Apass, Fs);
Hd = design(h, 'cheby1', 'FilterStructure', 'df2sos');

"Hd" at the interactive prompt gives me:

 sosMatrix: [1 2 1 1 -1.97232454484577 0.978340933957857;1 2 1 1 -1.94676435093801 0.948444394718458]
 ScaleValues: [0.00150409727802169;0.000420010945111697;0.891250938133746]

I notice that freqz(Hd.sosMatrix) gives me a wildly off magnitude estimation (but with the correct frequency cutoff).

My question is, how do I use this Hd.ScaleValues factor? There are 3 values for two sections.


texts: Which Approach Is Better for Decomposing an Image into High Frequency and Low Frequency Components?

texts:
Which approach is better or there is mathematical justification for using Bilater filter and Fourier Transform to decompose a image into High Frequency and Low Frequency Component.

Both Bilateral Filter and Fourier Transform can be used for getting High and Low Frequency component of an image.

Is one approach is better than the other or are both same in getting high and low frequency component of a image, or are they used in some particular context of image.

I know that Bilateral Filter uses range weights and Spatial weights into consideration for getting high and low frequency components.


texts: Best Way to Extract Frequency from Signal

texts:
I want to send signal that compound from 16 different frequency over GSM.

The problem is the GSM network is very noisy and I have only 400-3400 hz frequency available for use.

I try using fft algorithm for extracting the 16 frequencies but it can't handle the noise properly...

What is the best algorithm to go for my purpose?  I need:

being able to extract at least 16 different frequencies.
being able to do it fast (I must analyse at least 150 samples of 8
frequencies per second , and the goal is
do 450 samples of 16
frequencies per second

------- Added --------------------------

The goal is transport data over gsm (at 1200 bps minimum and the actual goal is 5200 bps), this mean that I need to modulate binary data over sound waves so it's will be transportable over gsm.

The idea is just use cellphone to transmit the data with dose encoded waves I created.

Gsm problems:

We have a lot of noise, we have 6 passes that each one add some noise (sound over trrs to cell phone, cell phone tocell tower , enter the gsm network, exit the gsm network, tower to cell phone , cell phone over trrs).
so the end results on the outer side is noisy.
The Telephony operators add filters on the gsm network to reduce noise, this filters looks on anything that dousn't look like human voice and trim it off the signal.
the bandwidth is 400hz-3400hz only so small bandwidth.

Now there is 3 main digital encoding methods over analog.
ASK, FSK, PSK.
AFK is not usable here because of the noise I described in the gsm section. (the amplitude disrupt severely)
PSK we have problems of sync and noise.
FSK looks like the most suitable for the this purpose because of it's ability to handle noise properly.

Now we have the issue of the bit-rate.
I want to use narrow fsk, so I will able to send at least a byte (8 bit) per sample (sample = the minimum samples we need to properly check which of the 8 frequencies exists in this part of the signal)
The problems I have is I cant demodulate properly the fsk over the GSM in a proper rate.

So the question is what is the best algorithm for a noisy environment to check the existence/absents of specific frequency in a signal where the importent parameters are:

the smallest amount of sample required from the signal (the smallest it is = I can send more byte's in a second)
the best suitable for the noisy environment.
can handle with the most frequencies (if I can send 32bit - 32 frequencies with the algorithm it's better then an 8 bit.... + keep in mind that it's full duplex = we receive and transmit the same amount of bit (in 8 bit case we have 16 frequencies in the signal 8 belong to the sender and 8 to the receiver) from each direction so 32 bit actual mean  is 32*2 frequencies on the signal where each side need only 32 bit (the sender read only the receiver frequencies and the receiver read only the sender frequencies).

What I did so far:

I check fft algorithm - didn't got me far...
looked into MUSIC algorithm
looked into the Goertzel algorithm.

There is probably a better algorithm to deal with this problem, hope some one here can point me to the correct direction


texts: Difference between MFB bound and ML bound for ISI channels?

texts:
I am considering a straight forward simulation setup in matlab.

I am transmitting a random sequence $\mathbf{x}$ of $\{+1,-1\}$ over a two-tap ISI channel $\mathbf{h}=[h_0 \hspace{2mm} h_1]$. The channel output is the noisy  filtered vector $\mathbf{y}=\mathbf{H}\mathbf{x} + \mathbf{n}$.
At receiver, I am implementing ML detector via Viterbi Algorithm to compute $\mathbf{\hat{x}}$.

I want to compare BER vs SNR for my simulation with the Matched Filtered Bound and Maximum Likelihood Bound. So what is the difference between these two and how can I compute these?
ML bound uses notion of $d_{min}$ but I don't understand how to compute it for any $\mathbf{h}$.

Can someone guide me here.


texts:
I had to work with sound(denoising algorithms), when I read about an algorithm it said that we do DTSTFT(Discrete-time short-time fourier transform) we will create an overlap between each "chunk" of data, for example if I have $$s=[1,2,3,4,5,6]$$ I divide it into $$s^{(1)}=[1,2,3,4]\\s^{(2)}=[3,4,5,6]$$ and then I use FFT on $s^{(1)}$ and $s^{(2)}$(this for $50\%$ overlap)

My problem is how can I do IFFT for this after I touched $S^{(1)}$ and $S^{(2)}$?  Everywhere I looked said to do the overlap but had not say a thing about the Inverse


texts: Estimating noise level without 'quiet' segment

texts:
I have a signal which consists of an onset and a release period.

The onset, sounding like a thud, excites frequencies across the entire frequency range (0-22 kHz). It lasts for about 0.1 s.

The release, sounding like a ring, excites 2-5 specific narrow frequency bands. Their amplitudes decay exponentially, albeit at slightly different rates.

I want to find the end of these narrow peak tracks. The method I am using now tracks the max amplitude in the narrow frequency bands and takes a windowed derivative. When this derivative passes through a threshold value, the end of the peak track is marked. This doesn't seem to be working consistently for particularly noisy signals.

My second idea was to track the amplitude of the peak tracks and compare them to the noise level. As soon as they reach the noise level, mark it as the end.

How do I estimate the noise level if there is no non-signal segment? The segment  before the strike, is often very short or has some ramp in amplitude making it unsuitable as an average.

Edit: Here's the spectrogram:


texts:
I know there are lots of papers proposing algorithms for the problem of reconstructing a signal from modulus of its Fourier Transform (so-called Phase Retrieval Problem).

Also, recently it is studied for the case of sparse signals.

I want to know whether there is any paper or survey working on "reconstruction of signal from modulus of its DFT" or equivalently, reconstruction of signal from its circular auto-correlation ?

Or is there a way we can relate these two problems?

Edit: The prior information I have is :

Time-domain signal is positive-valued.
Time domain signal is zero in some points, i.e. signal has a support $P \subset \{0,1,2,\cdots,N-1\}$


texts: How to extract noise when the noise on three data sets are near identical?

texts:
I have a microphone array with three microphones placed in a small equilateral triangle. During recording, there is a lot of noise introduced by an ADC. However, I am not interested in reducing this noise before sampling - for learning, I want to remove this noise digitally as best as possible, given that the noise on each channel is near identical.

Given that their noise is near identical, I imagine it would be possible to correlate the signals in some way, even in the presence of a stronger signal with slight time delays. If possible, how can I remove this noise?


texts:
I'm studying book about digital filter by Richard Hamming. And there is exercise to get odd and even expansion of g(x)=x where x is from 0 to $\pi$. I understood even expansion, but can't get into odd expansion:
$\pi$ - $2$*(${\rm sin} x$ +  $\frac{1}{2}$${\rm sin} 2x$ + $\frac{1}{3}$${\rm sin} 3x$...)

Can anybody explain odd expansion of this function g(x)=x?

Thanks in advance.

UPDATE:
Let me be more clear. The exersise is: show that g(x) = x has two expansions (0 < x < $\pi$)

$x$ =$\begin{cases}\pi - 2*(\sin(x)+\sin(2x)/2+\sin(3x)/3+...)
\\\pi/2 - 4/\pi*(\cos(x)+\cos(3x)/3^2+\cos(5x)/5^2+...)\end{cases}$

I merely understand cosine expansion, but I don't understand sine expansion...


texts: While converting discrete filter specs into continuous filter specs, why do we multiply gain by T

texts:
As I know so far that when we trying to design a DF(discrete filter) with specific specs "in terms of $\omega$", we take these specs and map it into continuous domain $H(j\Omega)$ using a certain method as "Impulse invariance" and design it by continuous design methods like "Butterworth".

As impulse invariance method says:
$$\omega = \Omega.T$$
and that part is fine to map $\omega_p(passband),\omega_s(stopband)$

But do we have to map also the gain of the filter?

i.e. $\delta_1,\delta_2$ and multiply them by $T$


texts: Will the capacity of a channel becomes unbounded if i increase its signal-to-noise ratio $S/N$ without limit?

texts:
According to the Shannon-Hartley theorem the capacity $C$ of a channel which has a signal-to-noise ratio of $S/N$ and a bandwidth $B$ can be calculated to be $C = B \log_2 \left( 1 + \frac{S}{N}\right)$.
Heren if $B→\infty$, the capacity does not become infinite since,
with an increase in bandwidth,the noise power also increases. If the noise
power spectral density is $N_0/2$,then the total noise  powr is $N=N_0B$,so the Shannon-Hartley law becomes
\begin{align}
C &= B \; \log_2 \left( 1 + \frac{S}{N_0B}\right) \\
&=\frac{S}{N_0}\left(\frac{N_oB}{S}\right)\log_2\left( 1 + \frac{S}{N_0B}\right) \\
&=\frac{S}{N_0}\log_2\left( 1 + \frac{S}{N_0B}\right)^{(\frac{N_0B}{S})}.
\end{align}
Now
$$\lim_{x\to0}(1+x)^{1/x}=e$$
So now it becomes
$$C_∞=\lim_{B\to\infty} C=\frac{S}{N_0}\log_2e=1.44\frac{S}{N_0}$$
so here channel capacity does not become infinite; that means, it's bounded.

But what happens if we increase the signal to noise ratio without bound? Will that give unbounded capacity? Is it possible?

After searching I got this but don't know how is this is possible; can some one justify this?


texts: Cross Correlation of output with input is equivalent to impulse response of a large bandwidth input (LTI System)?

texts:
What is the practical way to found the impulse response of a system ? As impulse function has a very large amplitude at zero time (analog signal).

So the cross corelation of output of a large bandwidth input signal (with relatively small amplitude) gives the impulse response of a LTI system .How is these two equivalent?


texts: Is "Windowed Fourier Transform" a synonym for "STFT"?

texts:
Sometimes I find the notation WFT (Windowed Fourier Transform) while other times I see STFT (Short Time Fourier Transform).
Are they the same?


texts: How to plot real time graph using python for the signals received from RTL SDR?

texts:
Using the below python code, I am able to plot a graph. Can anyone please tell   me how to modify this code to plot the graph continuously in real time?

from pylab import *
from rtlsdr import *

sdr = RtlSdr()

sdr.sample_rate = 2.4e6
sdr.center_freq = 93.5e6
sdr.gain = 50
samples = sdr.read_samples(256*1024)
sdr.close()

psd(samples.real, NFFT=1024, Fs=sdr.sample_rate/1e6, Fc=sdr.center_freq/1e6)
xlabel('Frequency (MHz)')
ylabel('Relative power (dB)')
show()


texts: How to separate Even and odd part (in-phase bit streams I and quadrature-bit streams Q) of the signal?

texts:
I am working on communication system in simulink and wondering, how can i separate the Even and odd part of the signal.
In my Transmitter block, I have :
$$
\boxed{\text{Binary Data signal}}{\longrightarrow}\boxed{\text{Serial to parallel}}^{\nearrow\boxed{\text{Even Part}}{\longrightarrow}\boxed{\text{Multiply PN code}}{\longrightarrow}\text{signal}}_{\searrow\boxed{\text{Odd Part}}{\longrightarrow}\boxed{\text{Multiply PN code}}{\longrightarrow}\text{signal}}
$$
Not sure how to separate those two part from signal ? Also can i use below equations for channels I and Q.
$$
I(t) = \sqrt{\frac{2E_s}{T} \cos(\omega t)\cos(45+n)}
$$
$$
Q(t) = \sqrt{\frac{2E_s}{T} \sin(\omega t) \sin(45+n)}
$$