texts: Failed to reproduce the smoothed spectrum plot from a time series

texts:
I tried to reproduce Watson's spectrum plot from these set of slides (PDF p. 30, p.29 of the slides), that came from this data of housing building permits.
Watson achieves a very smooth spectrum cruve in which it is very easy to tell the peak frequencies.
When I tried to run a FFT on the data, I get a really noisy spectrum curve and I wonder if there is an intermediate step that I am missing.
I ran the fourier analysis on python, using scipy package fftpack as follows:
from scipy import fftpack

fs = 1 / 12 # monthly
N = data.shape[0]
spectrum = fftpack.fft(data.PERMITNSA.values)
freqs = fftpack.fftfreq(len(spectrum)) #* fs

plt.plot(freqs[:N//2], 20 * np.log10(np.abs(spectrum[:N//2])))

Could anyone help me with the missing link?
The original data is:

Below is the Watson's spectrum curve, the one I tried to reproduce:

And these are my results:


texts: Interpreting the inverse fourier transform from a graph

texts:
I'm given a graph of the fourier transform of some function $x(t)$.  The graph is labelled $F(X(\frac{\omega}{\pi}))$ on the y-axis and $\frac{\omega}{\pi}$ on the x-axis.  The graph is plotted only on the interval from 0 to 1 (corresponding to $\omega$ values ranging from 0 to pi).  There are points at

$\left(\frac{\omega}{\pi},F(\frac{\omega}{\pi})\right)  = (0.2,1); (0.4,2); (0.6,3); (0.8,3);$ and  $(1,2)$

and is zero everywhere else.

I'm confident that these nonzero points correspond to a linear combination of unit sample delta functions centered at $\omega= \pi/5, 2\pi/5, 3\pi/5, 4\pi/5$, and $\pi$, respectively.  Since the delta functions are purely real valued, I conjecture that the inverse fourier transform of this function should be a linear combination of cosine functions.  Based upon a table of fourier transforms, I surmise that $$x(t)=\frac{1}{\pi}cos\left( \frac{\pi t}{5}\right) + \frac{2}{\pi}cos\left( \frac{2\pi t}{5}\right) + \frac{3}{\pi}cos\left( \frac{3\pi t}{5}\right) + \frac{3}{\pi}cos\left( \frac{4\pi t}{5}\right) + \frac{2}{\pi}cos\left(\pi t\right)$$.

Have I made the appropriate conclusions?


texts: What is the average of $\DeclareMathOperator{\rect}{rect} \rect(\cos(\pi t/2)) $?

texts:
We have this signal:
$$\operatorname{rect}\left(\cos\left(\frac{\pi t} {2}\right)\right) $$
I must find the average power , how can i get there ?
My solution:
I have seen that
$$-\frac 12 < \operatorname{rect}\left(\cos\left(\frac{\pi t}{2}\right)\right) < 1$$
is for
$$\frac 23+k \pi <t < \frac 43 + k\pi$$
so in a period it should be rectangular pulse , but how do I calculate the average ?


texts: Can spectral density be a complex quantity?

texts:
I have a signal ($S(t)$) which is product of a Gaussian ($G(t)$) and a random phase function ($e^{i\theta(t)}$, here $\theta(t)$ is a random function), as shown below

$S(t)=G(t).e^{i\theta(t)}$

If I calculate the auto correlation of such a signal ($E[S^*(t)S(t−τ)]$) it turns out to be a complex quantity and the same goes with the spectral density (Fourier transform of the auto-correlation function). My questions are the following.

If the above analysis valid? as the process is not wide sense stationary.
If the analysis is not valid, is there some way to handle this kind of situation?
If the analysis is valid, What does the complex spectral density signifies?


texts: I have seen that the SNR is measured for optical fibers "OSNR". What is the difference? Is there any equivalent?


texts:
Suppose I use the ideal lowpass filter eq to attempt creating a IIR filter using it
:
$$
h[n] = \frac{\sin(\omega_c \, n)}{\pi \, n}, \qquad -\infty < n < \infty
$$

Obviously, I can't really (unless I'm wrong) realize a filter with this impulse response because I don't have the finite difference equation to make the recursion happen.

But I'm told that it's also an unstable system (textbook solutions).
My question is.. why?

The stability condition states that a system is stable if its impulse response satisfies:
$$
\sum_{n=-\infty}^{\infty} \Big|h[n] \Big| < \infty
$$

From what I see, if I look at the condition blindly, it does satisfy the stability criterion. Doesn't it? Then why it's unstable? Is it because of something like a discontinuity at $n=0$ (even if technically its equal to one)?

Thanks!


texts: Correct way to build a basic filter bank (ie. for simple multiband audio processing)?

texts:
I need to build a filter bank for multiband processing. Last time I tried I was just starting at coding/DSP so I failed and didn't get a flat spectrum out of it. I just want to know the basic principle before I try again.

My understanding is you just need a series of high/low pass filters. The steeper the filters the more band separation you'll get.

I am guessing I would want a 12 dB/oct 2nd order design most commonly. So if I take some standard 2nd order high/low pass filters, like State Variable Filters, then I would want to set the Q to 1/sqrt(2) or 0.707 right? Because I don't want resonance.

Then I am confused as to whether I actually even need both high and low pass filters or if I'm just supposed to pick one type and build the whole thing with just that one type. ie. All low pass filters or all high pass filters.

To illustrate, let's say you wanted to split some white noise into 5 bands. You would have:

double whiteNoise = whiteNoiseGenerator();
SVFilter LPF1;
SVFilter LPF2;
SVFilter LPF3;
SVFilter LPF4;
SVFilter HPF1;
SVFilter HPF2;
SVFilter HPF3;
SVFilter HPF4;
float cutoffFreq1 = 250.f;
float cutoffFreq2 = 1000.f;
float cutoffFreq3 = 3000.f;
float cutoffFreq4 = 8000.f;
float filterQ = 1/sqrt(2);

Then you would set your filters in terms of Q and cutoff:

LPF1.setFilter(filterQ, cutoffFreq1);
HPF1.setFilter(filterQ, cutoffFreq1);
LPF1.setFilter(filterQ, cutoffFreq2);
HPF1.setFilter(filterQ, cutoffFreq2);
LPF1.setFilter(filterQ, cutoffFreq3);
HPF1.setFilter(filterQ, cutoffFreq3);
LPF1.setFilter(filterQ, cutoffFreq4);
HPF1.setFilter(filterQ, cutoffFreq4);

Then for processing what do you do next? Again, do you even need the HPFs? Couldn't you just subtract the prior filtered result from the next one? Like:

whiteNoise = whiteNoiseGenerator();
double outputBand1 = LPF1.process(whiteNoise);
whiteNoise = whiteNoise - outputBand1;
double outputBand2 = LPF2.process(whiteNoise);
whiteNoise = whiteNoise - outputBand2;
double outputBand3 = LPF3.process(whiteNoise);
whiteNoise = whiteNoise - outputBand3;
double outputBand4 = LPF4.process(whiteNoise);
whiteNoise = whiteNoise - outputBand4;
double outputBand5 = whiteNoise;

So you start with a full spectrum of white noise, and at each interval you just low pass to take a slice of that noise out and put it into the band you want. Since you subtract it from the remaining noise at each stage, the next stage only sees what remains from the prior filtrations. You don't need any HPF and the total amount of noise is unchanged since it all adds up so it must be flat.

I think that will work but it will create "slanted" crossovers between the bands along the axis of the LPF curves. To get straight crossovers I would need to use the HPFs too. Maybe like this:

whiteNoise = whiteNoiseGenerator();
double HPF1output = HPF1.process(whiteNoise);
double outputBand1 = LPF1.process(whiteNoise) - HPF1output;
whiteNoise = whiteNoise - outputBand1;
double HPF2output = HPF2.process(whiteNoise);
double outputBand2 = LPF2.process(whiteNoise) - HPF2output;
whiteNoise = whiteNoise - outputBand2;
double HPF3Output = HPF3.process(whiteNoise);
double outputBand3 = LPF3.process(whiteNoise) - HPF3output;
whiteNoise = whiteNoise - outputBand3;
double HPF4Output = HPF4.process(whiteNoise);
double outputBand4 = LPF4.process(whiteNoise) - HPF4output;
whiteNoise = whiteNoise - outputBand4;
double outputBand5 = whiteNoise;

Does that make sense? If so, do I set the cutoffs the same for each pair of high/low passes or do I need to stagger them somehow to keep things flat if it's say a 12 dB/oct 2nd order state variable filter?

Thanks


texts: Formant Analysis Frequency - What makes a peak value?

texts: from spectrograms, how can you tell which peak is really a formant. Like why is the pink arrow the first formant. And in the spectrogram below, which has fewer samples, what makes a value the peak value?


texts: How can I compute the inverse digital of a nonrational transfer function?

texts:
My $G(s)=1-e^{-s/\tau}$, $\tau$ is very small, say of order $10^{-4}$.

I need to compute a $H(z)$ (a digital filter) such that $H(z)$ has the inverse response of $G(s)$. Is ok even if $G(z)$ has the inverse response just for $f$<20k.

I don't not also zero and pole go $G(s)$.

I'm able to have an proximation of $H(z)$, but I found it with an iterative algorithm set $H(z)=\dfrac{\alpha}{1 -\beta z^{-1}}$.

I prefer a formal solution.


texts:
I am facing a trouble implementing the fft and appreciate your advice. My data structured as follow:
The data are retrieved from the share market for one day change in the share price. The time range from (11:00:00 to 15:00:00) this means 16200 Sec. In matlab its represented as serial so it will be from 0.4584 to 0.6458. Data are not sampled in fixed freqency. Price is sampled only when there is a change in price. For example ( it can be sampled in each second and sometimes in 30 Seconds interval).


texts:
I am trying to sample a sine wave and plot it's frequency components, but I am having problems implementing it.The result of taking 65536 samples of one cycle of a sine wave with max amplitude 1 and a frequency 100 can be seen below.Where the Y-axis this the magnitude of the complex Fourier sum, and the x-axis is the sample number.How can I see what frequency the sine wave has from this plot?, I would like to have a plot of freq vs Magnitude and not sample Number vs magnitude.

I was expecting one big spike of amplitude 1 or 2 spikes of amplitude 0.5 each, but I seem to be getting a large amplitude, I don't know what i might be doing wrong.

Matlab Code:

Amp = 1;
freq = 100;
dt = 2 * pi /65536;
index = 1;
for t = 0:dt:2*pi
   sine(index) = Amp * sin(2*pi*freq*t);
   sampleNumber(index) = index;
   index = index + 1;
end

transform = fft(sine);
magTransform = abs(transform);
plot(sampleNumber,magTransform);


texts: Wavelet signal analysis for ultrasound echo

texts:
I am new to the wavelet field, and I am working with ultrasound signal. For my study, I needed to measure material using the pulse echo reflection method.

Basically, when the transducer transmits the signal, the signal will reflected back to the transducer. If material has a thickness, the signal will reflected for two sides which are the front and back. After receiving both signals, I need to differentiate both of them, however the received signals are overlayed.

I think by using a wavelet method, I can differentiate the signal, but the result that I got is unclear.  Can anyone explain about this result?  Is it that both signals can be differentiated?


texts: What's wrong with this Whittaker-Shannon-Kotel’nikov interpolation implementation?

texts:
I tried to implement Whittaker-Shannon-Kotel’nikov interpolation formula but I get unexpected results: the reconstructed signal lags with respect to the original.

I know that I can not expect a perfect reconstruction because I am sampling a finite time signal (just one cycle of a sine wave) and this finite time signal is not bandwidth limited and so the Whittaker-Shannon-Kotel’nikov theorem does not hold... but I do not know whether to expect a lag.

I suspect there is something wrong in my implementation (for example in a previous implementation I was erroneously using $\frac{\sin x}{x}$ instead of $\frac{\sin \pi x}{\pi x}$) but I can not find it.

f=1;
T=1/f;
mysignal = @(t) sin(2*pi*f*t);

f_s=5*f; % sampling rate
T_s=1/f_s; % sampling period
t=0:T_s:T; % times at which mysignal is sampled
t1=0:.05*T_s:1.5*T; % times to plot mysignal "continuosly"

% plot mysgnal as a "continuos" signal
hold off;
plot(t1,mysignal(t1),'r')
hold on;

function y = Whittaker_Shannon_interpolation(samples,T_s,t)
% samples are the samples of the original signal;
% T_s is sampling period;
% t are the times at which interpolate the original signal;
% y are the interpolated values (reconstructed) of the original signal at times t.
  y = zeros(1,length(t));
  normalized_sinc = @(x) sin(pi*x)./(pi*x);
  for k = 1:1:length(samples)
    % https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula
    % the above formula uses so called "normalized sinc", i.e. sin(pi*x)/(pi*x)
    y+=samples(k)*normalized_sinc((t-k*T_s)/T_s);
  endfor
endfunction

samples=mysignal(t);
stem(t,samples,'b')
hold on;

y=Whittaker_Shannon_interpolation(samples,T_s,t1);
plot(t1,y,'k-o')

legend('original signal','samples','reconstructed')


texts:
do any of you know how to convert a FM index to a phase mod index equivalent considering they both have the same carrier and modulator frequency? Is it even possible?

This is in respect for sound synthesis and I'm talking about one sinusoid modulating another sinusoidal oscillator.

I'm being able to get very similar sounds with these parameters:

FM: carrier freq 200Hz, modulator freq 480Hz, index 400*

this means the carrier freq oscillates 400Hz up then back to original frequency and then down 400Hz and back 480 times a second

Phase Mod: carrier freq 200Hz, modulator freq 480Hz, index pi/4 **

** this index means the phase goes up 1/8 of a cycle and then back, down the same amount and back to the original point

FM and Phase Mod are clearly two independent processes, but I've seen people mention it can be equivalent, but I've never read anything that would tell me how to get the same results with both. So I don't know if it can really be done even if it sounds practically the same thing or if I just got lucky.

thanks


texts: Initial conditions for systems described in state space - LTI or not?

texts:
Suppose we have some system given by

$$\begin{aligned}
\dot{x}(t) &= Ax(t) +Bu(t) \\
y(t) &= Cx(t)+Du(t)
\end{aligned}$$

where $x(t)$ are the state variables, $y(t)$ is the output and $u(t)$ is the input. All matrices are constant. The same question applies to the discrete case

$$\begin{aligned}
x[n+1] &= Ax[n] +Bu[n] \\
y[n] &= Cx[n]+Du[n]
\end{aligned}$$

It is known that a system with non-zero initial conditions cannot be LTI. However, if $x(0)\neq0$, I don't see why the system above would not be LTI as it is expressed. As far as I know, if a system is expressed that way, it has to be linear and, as matrices do not depend on $t$, it should be time-invariant too.

So we have a system that has to be LTI as it is expressed in state space with constant matrices, but it can't be LTI because it has $x(0)\neq0$.

I can't see the mistake in the reasoning that leads me to this absurd contradiction. Can someone point it out?


texts: DC component in accelerometer data - filter before or after XYZ Euclidean sum

texts:
I'm working on filtering accelerometer XYZ data. The aim is gait-pattern recognition via finding characteristic frequency points and feeding these to a artificial neural network. The trouble is, and commonly seems to be so, filtering out a massive DC component (a strong peak at 0 Hz) in the power-vs-frequency domain, which overshadows higher frequency components. I'm working with Octave.

In literature, the typical approach is to gather accelerometer XYZ data, smooth some noise jitter with a sliding-window-average, apply a 4th order Butterworth highpass (at 0.9 Hz cutoff) and then take an Euclidean sum of the individual processed components (XYZ that is).

So far, whatever I try, there is always a strong peak centered at 0 Hz. Trying has included:
1) Subtracting mean of raw data (x=x-mean(x), for all channels)
2) Various orders of Butterworth on individual channels before the Euclidean sum, in various settings

What did get rid of the strong 0 Hz peak was applying the Butterworth filter on the Euclidean sum of XYZ channels. The characteristic 2 Hz gait peak remains and some higher order frequencies are amplified.

My question is: is there anything fundamentally wrong in filtering the XYZ data after the Euclidean sum? Or: is this correct? Or: what is happening, why is the 0 Hz peak not getting filtered with individual-channel filtering (that would be a hint at something non-distributive).

Here are some plots:
Raw XYZ accelerometer data
Raw XYZ accelerometer data http://www.hot.ee/jaaniussikesed/rawdataxyz.png

Data with very slight sliding window average (pretty much raw)
Data with very slight sliding window average (pretty much raw) http://www.hot.ee/jaaniussikesed/walk_xyz_nf_fft.png

4th order Butterworth, cutoff at about 1 Hz, on individual XYZ before Euclidean sum
4th order Butterworth, cutoff at about 1 Hz, on individual XYZ before Euclidean sum http://www.hot.ee/jaaniussikesed/walk_xyz_f_fft.png

4th order Butterworth, cutoff at about 1 Hz, taken on a Euclidean sum of XYZ accelerometer data. Notice no 0 Hz peak.
4th order Butterworth, cutoff at about 1 Hz, taken on data after Euclidean sum of XYZ http://www.hot.ee/jaaniussikesed/walk_euc_f_fft.png

4th order Butterworth filter frequency response:
4th order Butterworth filter frequency response http://www.hot.ee/jaaniussikesed/frq_response.png

Some data-data:
Samples: 1112
Sample rate: 84 samp/s

Butterworth implementation in Octave:

Fc= 0.5; #Cutoff frequency in Hz
Wcn=Fc/(srate/2); #Normalized cutoff frequency
order=4; #Order of filter
[B,A]=butter(order,Wcn,"high"); #Define Butterworth highpass filter

Note1: changing the cutoff frequency between 0.1-1 does not really change the picture.
Note2: there have been quite many posts on a very similar problem, but they all have their nuances, as does this I believe.

Thank you for any help or/and advice!

EDIT1:
Here is the code (running Octave 4.0.0, packages: control, data-smoothing, optim, signal, struct):

# Load data
load("ID:walk-2016-05-02T22:58:23+02:00.mat");

# Separate XYZ and clip
clipc=300;
ax=datamat(clipc:(end-clipc),1);
ay=datamat(clipc:(end-clipc),2);
az=datamat(clipc:(end-clipc),3);

srate=floor(srate); #sample rate

# Design filter
pkg load all
Fc= 1; #Cutoff frequency in Hz
Wcn=Fc/(srate/2); #Normalized cutoff frequency, srate (sample rate) loaded from file
order=4; #Order of filter
[B,A]=butter(order,Wcn,"high");

# Apply filter
axf=filtfilt(B,A,ax);
ayf=filtfilt(B,A,ay);
azf=filtfilt(B,A,az);

#Euclidean norm
afen=(axf.^2+ayf.^2+azf.^2).^(1/2);

#FFT
N_fft=2^13; #Frequency sample points, empirical
afft=fft(afen,N_fft); #Centered double sided FFT
freqpwr=afft.*conj(afft)/(N_fft*rows(afft)); #Power of each frequency spectrum component
freqpoints=srate*(0:N_fft/2-1)/N_fft; #Bins/points of FFT frequency axis

#Plot power spectral density
plot(freqpoints,freqpwr(1:N_fft/2));
title("One Sided Power Spectral Density");
xlabel("Frequency (Hz)");
ylabel("Power spectral density");

fprintf("Program paused. Press enter to continue.\n");
pause;
close;

Link to data file.


texts:
I am new to this site and field too.I always thought analog is continuous and digital signal is discrete. I read this today and got confused about analog signal

  Analog can be
  continuous time(CT) or discrete time(DT). A discrete time signal is not a digital signal; a
  discrete time signal is one in which the amplitude is a continuum, and not discretized. On
  the other hand, if the time is discretized, this signal is still analog. Only when a discrete
  time signal is passed through an A to D converter, it becomes a digital signal.

Will anybody explain me what actually it means?


texts: Effects of linear interpolation of a time series on its frequency spectrum

texts:
Situation
In order to synchonisize different time series i have to apply linear interpolation on them. After the interpolation and synchronization the
signal is transferred into its frequency domain for further analysis.
The interpolation step should therefore not alter the frequency domain information.
Experiment
Via python I created a signal of white noise.
The second signal is based on the first but shifted by half of the original sampling interval. The third signal was created by oversampling the original signal by the factor 5.0.
Afterward all three signals were transferred to the frequency domain (Bottom of the Figure)
Questions

Why does the curve of the shifted signal decline for higher frequencies?
Why does the curve of the over sampled signal decline for higher frequencies? Though the over sampled signal has the exact same shape as the original signal.
What can i improve in order to interpolate the original signal at specific points but minimize the effects on its frequency spectrum?

My python code for the experiment and the plot: https://pastebin.com/iz08Hiud
Many thanks in advance!
Update
Thank you all for your comprehensive adn helpfull answers, especially robert bristow-johnson and howpow2 for directing me on the right path. I am using now an implementation of the sinc kernel interpolation, which shows an almost perfect frequency response.
See the bottom row of the plot for the results in form of the error between the original frequency spectrum and the frequency spectrums of the interpolated signals


texts: I have a Gaussian random variable $z$, whose mean $\mu$ and variance $\sigma^2$ are functions of two variables $x$ and $y$. Is there any way in MATLAB so that I can find the values of $x$ and $y$ for which $P(z)$ = $N(z$;$\mu$$(x,y)$,$\sigma$$(x,y))$ is maximised for a given value of $z$? Thanks in advance.


texts: Identify large solid color areas of image

texts:
This is a simple image classification problem. I want to categorize jigsaw puzzle images to difficulty levels (e.g, easy, medium, hard).

I want to develop an simple algorithm for this classification, which is based on factors like large solid colored areas, repeating patterns on the image (more details). Such classification can be seen in following website.
example 1

I want to automate this without much human interventions. Is this possible and whats the best way to approach it. (not worried about performance, simple is good, does not need huge accuracy)

EDIT

I was able to generate the magnitude and phase plots of the 2d FFT. However, I am not clear how to come up with a classification by looking at magnitude plot. Any ideas.


texts: Dimensional analysis in the Morlet wavelet

texts:
I'm a complete beginner on signal processing, so please forgive the basicness of this question - but I'm very confused by this. Regarding the Morlet wavelet, Wikipedia says:

  The wavelet is defined as a constant $\kappa_\sigma$ subtracted from a plane wave and then localised by a Gaussian window:
  $$ \Psi_\sigma(t) = c_\sigma\pi^{-\frac{1}{4}}e^{-\frac{1}{2}t^2}(e^{i\sigma t} - \kappa_\sigma) $$
  where $\kappa_\sigma=e^{-\frac{1}{2}\sigma^2}$ is defined by the admissibility criterion, and the normalisation constant $c_\sigma$ is:
  $$ c_\sigma = \left( 1 + e^{-\sigma^2} - 2e^{-\frac{3}{4}\sigma^2} \right)^{-\frac{1}{2}} $$
  ......... The parameter $\sigma$ in the Morlet wavelet allows trade between time and frequency resolutions. Conventionally, the restriction $\sigma>5$ is used to avoid problems with the Morlet wavelet at low $\sigma$ (high temporal resolution).

Presumably, given the comment about the typical restriction of $\sigma>5$, this trade-off parameter $\sigma$ is meant to be a dimensionless value. But then (or, in fact, whatever the dimension of $\sigma$ is), the formula for $\Psi_\sigma(t)$ is not dimensionally invariant under changes of the units in which time is measured! Shouldn't there be at least some extra parameter in front of $t$, so that the formula makes sense in terms of dimensional analysis, if it is meant to be applicable to real-world data (as suggested by comments about "conventional restriction" regarding trade-off)?

(Incidentally, I've seen elsewhere much simpler formulae for the Morlet wavelet, which do make dimensional sense, and where $\sigma$ has the dimension of time - which makes sense, as low $\sigma$ means a narrow time-window of significant relevance.)


texts:
I am practicing for an exam and I need help in the following question (it is just for me- no grade...)

we get an input image, 6 convolution filters and 6 results.
we need to find the matches between the filters to the result images.

so far I found that filter 3 matches to image D and filter 4 matches to image f.

now I stuck and I would appreciate any help.
thanks!

EDIT: the origin of all the filters is in the center.


texts: Power spectral density vs Energy spectral density

texts:
I read the following on Wikipedia:

  Power spectral density:

  The above definition of energy spectral density is most suitable for
  transients, i.e., pulse-like signals, for which the Fourier transforms
  of the signals exist. For continued signals that describe, for example, stationary physical processes, it makes more sense to define a power spectral density (PSD), which describes how the power of a signal or time series is distributed over the different frequencies, as in the simple example given previously.

I don't quite understand that paragraph. The first part says that "for some signals .. the Fourier transform does not exist".

For which signals (in the context that we are discussing) does the Fourier transform not exist, and we therefore need to resort to the  PSD  rather than using the energy spectral density?
When obtaining the power spectral density, why can't we compute it directly? Why do we need to estimate it?
Finally, on this topic, I have read about methods that use Kayser-windows when computing the PSD over time. What is the purpose of these windows in PSD estimation?


texts: Why Is the Canny Edge Detection Used Instead of Sobel / Prewitt Edge Detection Before Hough Transformation?

texts: I know that for the Hough Transform to work on an image, it needs to be a binary image. To convert from a grayscale image, an edge detection algorithm should be employed. I notice that people always use Canny edge detection instead of others (Sobel, etc). Why is that?


texts: Is there a way to apply blind source separation using ONE mixture?

texts:
This question is a follow-up for my previous one: Reduce the "glare" effect of a receiver within a system.

I have broken the problem into pieces and I humbly think that I've got a way to solve it it with Duet BSS. However, there is still a problem: I have got

$$Y(t)=h_1(t)\otimes X(t) + h_2(t)\otimes s(t)$$

Is there a way be it a hardware or a software solution to use ONE mixture to get my two signals of which I only know that $X$ is much stronger than $s$ and that on the spectrogram they don't overlap?

I thought about putting a relay mast between the emitters and the receptor so I can get two rich mixtures. However, there is a risk that I increase, by doing so, the "blinding" effect on the receptor.

The signals are not periodic but could be seen as cyclostationary.


texts:
We are analyzing a real signal with the DTFT. Since we are using a limited number of samples it's like we are transforming a finite signal.

As I remember, the FT of a finite signal has an infinite band spectrum. So I suspect that the DTFT of a sampled signal (with a finite number of samples) always produce aliasing (since in the DTFT spectrum each replica is an infinite band). Is this correct?


texts: Butterworth-like filter with no zero at Nyquist?

texts:
I have an old digital synthesizer that has a Butterworth-like 2nd order filter with resonance. However, with high frequency tones it appears there is no steep cutoff near Nyquist (I am sure it is the Nyquist as the tones fold beyond that point), as is usual with the bilinear-transformed Butterworth filter. This poses few questions.

Is the attenuation caused by aliasing or something different, and what is it called?
Are there digital alternatives to the Butterworth or bilinear transform that do not suffer from the attenuation near Nyquist, or is this a usual feature?
What filter could an old resource-strained synthesizer have used to achieve a response without the attenuation?

Clarification: By Butterworth-like response it is meant that the response appears similar at lower frequencies, I have not confirmed that there is no ripple or anything like that. Nonetheless, a 2nd order filter of some sort.


texts: Trying to understand WAVELET TRANSFORM Frequncy-Time diagram

texts:
I am reading

The Wavelet Tutorial

Part III
MULTIRESOLUTION ANALYSIS
&
THE CONTINUOUS WAVELET TRANSFORM

by
Robi Polikar

The author explains about the following fig and says;

  Note that boxes have a certain non-zero area, which implies that the
  value of a particular point in the time-frequency plane cannot be
  known. All the points in the time-frequency plane that falls into a
  box is represented by one value of the WT.

I am unable to interpret this fig and also that why value of a particular point in the time-frequency plane cannot be known.   I can always intersect two lines one passing through a particular value of t and another through f. The intersecting point will give me exact location of frequency at particular time (and V.V). I am not able to understand what thsi diagram is all about and what is its significance?


texts: Proof $GM$ (in dB) of $2nd$ order system is $\infty$

texts:
How can we prove that the Gain Margin $[GM]$ of $2nd$ order system is $\infty\quad ?$

My Approach:
Let us consider a $2nd$ order open loop system : $$G(s)= \frac{k}{(s+1)(s+2)}$$
Now, we know: $GM(in\quad dB)=-20log_{10} |G(jw_p)|$
$$ For \quad w_p \quad : \quad \angle GM(jw)=-180$$  $$\implies -tan^{-1}(w)\quad-tan^{-1}(w/2)=-180$$ $$\implies tan^{-1}(w)\quad+tan^{-1}(w/2)=180$$ $$\implies w=0 \quad or \quad w_p=0$$ $$\therefore |G(jw_p)|=k/(1 \times 2)=k/2$$ $$\implies GM(in\quad dB)=-20log_{10} (k/2) \neq \infty$$

so where is my mistake and how can i proof this results....??please help


texts: Why does my zoom audio, when there are nearby microphones, seem to usually take an open loop format?

texts:
The signal from the speaker on device A gets into the microphone from device B, then the amplifier and then the speaker from device B, back into device A. The result is a sort of runaway signal that is always reamplifying itself?
Why can't I create a sort of closed loop system when I have nearby microphones? What would I need to do? What gear and/or software would I need to create a closed feedback loop?


texts: Single-pole IIR low-pass filter - which is the correct formula for the decay coefficient?

texts:
A single-pole IIR low-pass filter can be defined in discrete time as y += a * (x - y), where y is the output sample, x is the input sample and a is the decay coefficient.

However, the definition of a varies.
On Wikipedia, it's defined as 2πfc/(2πfc+1) (where fc is the cutoff frequency).

But here, a is defined as follows: 1 - e^-2πfc.

Their graphs look similar, but which one is more accurate?

The blue graph is the Wikipedia formula, the green is the second one, fc is the x-axis.


texts: How to scale audio signal from PCM data to a dB SPL value?

texts:
I've got a set with MEMS microphones that measure audio signals at ultrasonic frequencies. The MEMS microphones have a PDM output which is then converted into PCM (this is necessary to allow for a microcontroller to do certain processing on the sampled audio data).
I'm trying to come up with a method to convert the PCM samples into dB SPL and the best resource I've found on this is this link: https://curiouser.cheshireeng.com/2015/01/16/pdm-in-a-tiny-cpu/. I understand how they calculate a RMS value from 977 PCM samples (this is called an SPL value in internal logarithmic units in the article). They relate this RMS value to a dB FS value by using the microphone datasheet (where the maximum possible PCM value/RMS value for a square wave will be equivalent to a known maximum value of dB FS of +3 dB FS ). I don't understand how the author then creates a linear relationship between dB FS and dB SPL (akin to the classic y=mx+b). The specific paragraph from the article discussing this is listed below:

To relate the finished SPL value in internal logarithmic units to dB SPL, we have to note that the microphone data sheet claims that a 1 kHz tone at 94 dB SPL will typically register as -26 dB FS, where 0 dB FS is largest amplitude sine wave that can be represented in a PCM sample without clipping. A full scale square wave is then +3 dB FS, and which would measure as 8 * log2(8192) or 104 which can be rescaled to dB by multiplying by 20 * log10(2) / 8 or about 0.75 to get 78. Subtract the 3 dB for a sine wave, and we find that the offset is about -75 dB to dB FS, or +19 to dB SPL.

Putting this all together, if we wanted to output true db SPL we would need the following expression in terms of our computed variable spl: dB SPL = (3 * spl / 4) + 19

I'm not understanding how the coefficient of 0.75 or the intercept of +19 is justified. Does anyone have any ideas or some additional resources I can consult on this?


texts: Fourier Transform for frequency analysis

texts: I applied Fourier transform on a sound (PCM generated digitally with sin function) and frequency detection is accurate. However, if I play the sound and record it, and then apply a Fourier transform to detect frequencies, my results are nowhere close to original results. The window size I have chosen for Fourier transform is 30 samples. ( Sampling frequency is 44100 Hz and 30 samples represent a symbol. Each symbol frequency is shifted by 1470 Hz.). Is my approach right ?


texts: How to calculate the mean/center frequency of the spectrum?

texts:
Assume that the Fourier Transform of $x(t)$ is $X(j\omega)$. And I want to calculate the $\bar{\omega}$ which is the center freqency in the spectrum ( the highest).

And a classmate tells me  this formula $$ \bar{\omega} = \int \omega X(j\omega)\ d\omega$$ could be used to calculate the center frequency ? But it is the first moment, it seems that the result is mean frequency.
Could anyone help me with this question?
What is the mean and center frequency of the spectrum?

There is no solution in Alan's book.

EDIT:
Maybe I did not state the center frequency clearly. I want to find out the frequency where the amplitude spectrum is biggest.
Another situation: assume that the signal is discrete $x[n]$. So we could just compute a period $[-\pi, \pi]$. In this case, this formula could represent the center frequency?


texts:
I'm using Math.net and Bass.Net to write FX effects for guitars. I'm currently stuck on the the pitch shifter. I have to implement a custom DSP callback function which provides the guitar PCM data and the length in bytes (the sample rate is 44100 Hz). I currently use the FFT implementation of Math.Net, but I don't know what to do with the complex output. My code is the following:

        //buffer (IntPtr) contains the PCM data in 32-bit floating point samples
        //length (int) is the number of bytes to process

        Complex[] cOutput = new Complex[length/4]; //allocate the buffer to store the complex numbers
        float* data = (float*)buffer;
        for (int i = 0; i < (length / 4); i++)
        {
            cOutput[i] = new Complex((double)data[i], 0.0); //fill the Complex array
        }

        Fourier.Forward(cOutput); //apply the FFT to the filled Complex array and overwrite its data

So I'm now have an array of Complex, but I don't know what part of this array I have to modify to change the pitch. So what comes next in the pitch shifting algorithm?

Thanks for your help


texts: Techniques for detecting small blobs in noisy images

texts:
I am trying to write a program that uses computer vision techniques to detect (and track) tiny blobs in a stream of very noisy images. My data is of two types: one set of images are not so noisy, which I am just using to try different techniques with, and the other set are noisier, and this is where the detection needs to work at the end. Here are some samples of the higher-noise images. The blobs that need to be detected are the small black spots near the center of the image.

Initially I started off with a simple contour/blob detection techniques in OpenCV, which were not very helpful. Eventually I moved on to techniques such as "opening" the image using morphological operators, and subsequently performing a Laplacian of Gaussian blob detection to detect areas of interest. This gave me better results for the low-noise versions of the images, but fails when it comes to the high-noise ones: gives me too many false positives. Here is a result from a low-noise image (please note input image was inverted).

The code for my current LoG based approach in MATLAB goes as below:

while ~isDone(videoReader)
    frame = step(videoReader);
    roi_frame = imcrop(frame, [660 410 120 110]);

    I_roi = rgb2gray(roi_frame);
    I_roi = imcomplement(I_roi);
    I_roi = wiener2(I_roi, [5 5]);
    background = imopen(I_roi,strel('disk',3));

    I2 = imadjust(I_roi - background);
    K = imgaussfilt(I2, 5);
    level = graythresh(K);
    bw = im2bw(I2);

    sigma = 3;
    % Filter image with LoG
    I = double(bw);
    h = fspecial('log',sigma*30,sigma);
    Ifilt = -imfilter(I,h);

    % Threshold for points of interest
    Ifilt(Ifilt < 0.001) = 0;
    % Dilate to obtain local maxima
    Idil = imdilate(Ifilt,strel('disk',50));

    % This is the final image
    P = (Ifilt == Idil) .* Ifilt;

Is there any way I can improve my current detection technique to make it work for images with a lot of background noise? Or are there techniques better suited for images like this?


texts: Distinguishing FIR and IIR from difference equation

texts:
Find the transfer function of the difference equation $$y_n = x_n + 1.2y_{n-1}$$

I fail to understand how one can distinguish between an FIR filter and an IIR filter by looking at the equation given above.

How to find if an IIR or a FIR filter is stable from the transfer function?

$$ H(z) = \frac{y(z)}{x(z)}$$

$$x_n = y_n-1.2y_{n-1}$$

$$x(z) = y(z)-1.2z^{-1}y(z)$$

$$H(z) = \frac{y(z)}{y(z) -1.2z^{-1}y(z)}$$

$$H(z) = \frac{z}{z-1.2}$$

Thus we arrive at the transfer function.

I am a week into signal processing. Any suggestions and tips appreciated. Please forgive my mistakes.


texts: Find integral of DTFT after sampling (Graph of CTFT given)

texts:
So for the first question:

If this is sampled at 10kHz, then the amplitude is scaled by 10000. In the DTFT, the frequency 3.5kHz gets mapped to 3.5/10* 2pi=0.7pi. So this point lies outside the range of the integration.

The next non-zero point is, 10-3.5=6.5kHz. This gets mapped to 6.5/10*2pi=1.3pi. This also lies outside of the integration range.

So the integral in the range 0.75pi to pi contains only zero points, which means the integral is zero. Is this correct?

For the second question,

The frequency 3kHz gets mapped to 3/8*2pi=0.75pi. This is the lower limit of the integral. 3.5kHz gets mapped to 3.5/8*2pi=0.875pi. So the spike at 3.5kHz is inside the integral. The next non zero point, 8-3.5=4.5kHz gets mapped to 1.125pi, which is outside the integral.

So in this case, the integral from 0.75pi to pi contains only one non-zero spike at 0.875pi of amplitude 8000. However, I don't know what is the integral of a spike. Please help.


texts: I am experimenting with some signals. I ma generating a 2.47 Ghz signal from signal generator and receiving it by an SDR (LMS7002m Rf chip) with Rx frequency tuned to 2.47 GHz. Upon plotting the received I/Q data I get the following plot. Can anyone help me understanding why there is a second circle below the circle centered at (0,0). Ideally I should be seeing a single circle..


texts: Can the instantaneous frequency be always derived from an analytic signal?

texts:
Instantaneous frequency can be defined as a derivative of an instantaneous phase of an analytic signal which can be nicely seen in practice in this example from Scipy's documentation. But it seems like it does not always work like this. I played with the code from the example and obtained varied results like for this sine wave with frequency gradually changing from 2 to 12Hz:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert, chirp

duration = 1.0
fs = 400.0
samples = int(fs*duration)
t = np.arange(samples) / fs
w = 2*np.pi*(2 + 10*t)
signal = np.sin(w*t)
analytic_signal = hilbert(signal)
amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase) /
                           (2.0*np.pi) * fs)

This is the resulting figure with the original frequency plotted in grey on the bottom subplot:

fig = plt.figure()
ax0 = fig.add_subplot(211)
ax0.plot(t, signal, label='signal')
ax0.plot(t, amplitude_envelope, label='envelope')
ax0.set_xlabel("time in seconds")
ax0.legend()
ax1 = fig.add_subplot(212)
ax1.plot(t, 0.5*w/np.pi, lw=2, color='gray')
ax1.plot(t[1:], instantaneous_frequency)
ax1.set_xlabel("time in seconds")
fig.show()

In the bottom plot the blue line and the grey line go in the same direction, but there is a huge discrepancy between the two and I do not mean the ever-present edge effect. The calculated instantaneous frequency is rising at the rate that is almost a double of the actual reaching values above 20Hz at the end.

I tried it for other sine waves with constant amplitudes and varied frequencies getting similar results. Only when I used a signal with a constant frequency the blue line was matching the grey line on the bottom plot. My question here is: are there any conditions a signal must meet for the definition of the instantaneous frequency to hold true, so it can be derived from the instantaneous phase of the analytic signal?


texts: What is the definition of the bandwidth of a signal?

texts:
I was reading over the famous Nyquist–Shannon sampling theorem and I was trying to understand what the term bandwidth meant rigorously (exactly) with no ambiguity. In this site the define it as follow:

  2) Bandwidth is the range of frequencies -- the difference between the highest-frequency signal component and the lowest-frequency signal component -- an electronic signal uses on a given transmission medium. Like the frequency of a signal, bandwidth is measured in hertz (cycles per second). This is the original meaning of bandwidth, although it is now used primarily in discussions about cellular networks and the spectrum of frequencies that operators license from various governments for use in mobile services.

which to me it means that given some periodic function:

$$ f(x) = \frac{1}{2}a_0 + \sum^{\infty}_{n=1} a_n cos(nx) + \sum^{\infty}_{n=1} b_n sin(nx) $$

we just need to extract the frequencies that are the largest and subtract them. For the sake of an example consider:

$$ f(x) = cos(2 \pi x) + cos(4 \pi x ) + sin( 3 \pi x) + sin( 7 \pi x) $$

the frequency of a sin/cos is $f = \frac{n}{2 \pi}$, thus, components of $f$ with the largest and smallest $n$ give the frequencies we want. So the band width $B$ of $f$ is:

$$ B = f_{max} - f_{min} $$

$$ B = \frac{n_{max} }{2 \pi} - \frac{n_{min}}{2 \pi} = \frac{7 }{2 \pi} - \frac{2}{2 \pi} $$

is this algorithm correct? Do I have the correct understanding of bandwitdh in the context of signals and systems?

I think less like an engineer and more like a mathematician so it was important for me to know exactly what bandwidth means in the context of Fourier series since its a totally unambiguous example to me (since all period functions can be expressed like that).

Also, as I was about to tag the question the definition of band width according to this site came up as:

  the difference between the upper and lower frequencies in a contiguous set of frequencies.

this reminded me that a lot of the time I see signals represented in what this community calls "frequency domain", as sometimes a continuous line (i.e. $f(c) = lim_{x \rightarrow c}f(x)$). This seems to be really unintuitive to me because the Fourier series is a countable summation but a continuous frequency domain plot would imply its a uncountable summation of the periodic function. In the light of this, how is "contiguous set of frequencies" a good definition? Or what does "contiguous" mean?

As you can see I just want to understand what bandwidth means in the end in a mathematical sense.


texts: Can thresholding in wavelet transform be implemented on the approximation coefficients?

texts:
Sorry as it might be a silly question since I am a beginner in wavelets. I have been using the matlab wavelet toolbox for denoising purposes and saw that thresholding is only applied on the detail (high pass) coefficients. Is there a different method that can be used where thresholding is applied on the approximation coefficients rather than the detail coefficients ?
I am interested in this because I want to implement the denoising in real time and would like the least amount of levels for decomposition for low computational power, and most of my noises fall in the low frequency range so thresholding on approximation coefficients would be helpful.


texts: How we can obtain the time-integration of a sine signal?

texts:
could we calculate the below integral by the Fourier series or the Fourier transform properties?

$$\int_{-\infty}^t \sin(\omega_0\tau)d\tau=?$$


texts: Some software (for example SDR Sharp) has a "Correct IQ" option.  This is used to get rid of the DC spike in the middle of the display.  What does this setting do to the signal?  What is wrong with the IQ values that needs to be corrected?


texts: amplitude of upsampled and downsampled signal without filter

texts:
Given:
$$
DTFT\{x[n]\}=X(\omega)= \begin{cases}
1 & |\omega| \leq 2/\pi
\\
0 & 2/ \pi < |\omega| < \pi
\end{cases}\ \ \ \ \ (periodic\ 2\pi)
$$

If I downsample $X(\omega)$ by M. I get:

$$
\begin{aligned}
x_d[n]      &= x[nM] \\
\\
X_d(\omega) &= 1/M\ X(\omega/M)
\end{aligned}
$$

Now what happens if I up-sample $x[n]$ again by $L$?

why can't I just let $M = 1/L$ and substitute into down-sampling equation to get up-sampling equation?

$$
\begin{aligned}
x_i[n]      &= x[n/L] \\
\\
X_i(\omega) &= L\ X(L \omega)
\end{aligned}
$$

I'm a little bit confused because my book is telling me that upsampling equation doesn't scale the amplitude by L:

$$
X_i(\omega) = X(L\omega)
$$

Just wondering what happened to the L term when upsampling?

(page 112, Schaum's Outline, Digital Signal Processing, Second Edition)

If I look up DTFT time scaling property, I find:

$$
DTFT\{\ x[an]\ \} = 1/a X(\omega / a)
$$

Which seems to confirm that upsampler should have an L term ..


texts: Why does FFT output need to conjugate processing?

texts:
I came across some matlab code as the below,

fftout   = fft(adc_in, 1024);
fftout_p = fftout.*conj(fftout);

I think fftout is formatted as complex form. from here I'm ok.

Next stage that signal is multiplicated by itself which is conjugated.

I think it means that it want to get the real number not complex number.

From here, I'm not sure and want to know that practically what does Imaginary part affect to real signal processing?

Can we ignore like above that?

In real practically engineer field, What do conjugating and Imaginary part mean and affect to the real engineer field?


texts:
I have a signal, shown in the attachment (blue one). As you can see that the latter part of my signal has a higher dc offset. I want to bring the signal into the same level. I have tried eliminate the low frequency signal by using FFT (I FFT the signal, then set the low frequency element to be 0, then returned it by IFFT), but the result that I get is shown in red color.
I am open to any suggestion and any python code is very appreciated.

Additional info:
sampling frequency = 100 Hz
Low cutoff frequency = 0.05 Hz


texts: Minimum Sampling frequency for this signal to avoid aliasing

texts:
I would like to compute the minimum sampling frequency for the signal $x(t)$ given below:

$x(t)= sinc ^3(300t) * sinc ^5(600t)$ and  $ *$ stands for convolution.

where the sinc function is defined as $sinc (w t) = sin (\pi wt)/(\pi wt)$.

I have tried like this  : $sinc^3(300t)$ will correspond convolution of 3 rectangular functions in frequency domain which results in a frequency limited function to 3 times 150 Hz ie 450 Hz. And the other function corresponding to maximum of frequency 5 times 300 Hz which is 1500 Hz. Now finally convolution operator get transformed to multiplication in frequency domain , is the maximum frequency content limited to 450 Hz and minimum sampling frequency 900 Hz?

It would be helpful if someone shows me the derivation of fourier transform of sinc function for this purpose.


texts: Minimum Number of Errors required for BER Curves in Matlab Simulation

texts: what should be minimum Number of Errors required for getting BER Curves in Matlab Simulation for different number of transmitted bits ?


texts: Why do I have to use variance and not mean value to calculate noise?

texts: Why, when I have to calculate some types of noise (for example shot noise and others which are in the denominator of SNR ratio) or attenuation from turbulence, I have to use variance of the of the density probability and not the mean value? In all of the formulas they use  variance (the second moment) and this is not clear for me.


texts: Discrete wavelet transform disadvantages

texts:
I read in a paper that the  discrete wavelet transform (DWT) has two disadvantages
The first one is the shift variance property due to the downsampling process. Could you please help me understanding why downsampling leads to shift variance?
The second disadvantage is the lack of directional selectivity.  Why does DWT lacks directional selectivity?


texts: Separate speech signals that are bandlimitted

texts: I have 2 kind of signals (refer figures) in a folder. One is band limited to 4khz and another is band limited to 8khz, both are with sampling frequency 16khz. I need to separate them. Any idea??


texts:
I asked this at Stack Overflow, they told me to come here.

First, please excuse my ignorance in this field, I'm a programmer by trade but have been stuck in a situation a little beyond my expertise (in math and signals processing).

I have a Matlab script that I need to port to a C++ program (without compiling the matlab code into a DLL). It uses the hilbert() function with one argument. I'm trying to find a way to implement the same thing in C++ (i.e. have a function that also takes only one argument, and returns the same values).

I have read up on ways of using FFT and IFFT to build it, but can't seem to get anything as simple as the Matlab version. The main thing is that I need it to work on a 128*2000 matrix, and nothing I've found in my search has showed me how to do that.

I would be OK with either a complex value returned, or just the absolute value. The simpler it is to integrate into the code, the better.

Thank you.


texts: How to choose the number of codewords in bag of words model

texts: In the bag of words model, generally is it better to have more codewords or less? For example, if I have 1000 x 20 raw features, then what would be good value for codebooks? I am thinking of smaller size codebook ( such as 20) so that it forces the representation to be sparse.


texts: Can I use a filter design to optimize the preconditioning of an optimizer?

texts:
Consider a noisy time series y_i and that I have waited until end of experiment. I now want to fit a nonlinear parameterized function F(A,d,k; t) to determine A,d,k.

So I can happily Newton-Raphson my way through | y - F|^2. But gosh the y's are noisy. I can't stop myself from filtering them first to throw away some high frequency junk. I think it helps. So I am minimizing |L y - F|^2 ... All is good. Except. To design L, I did some handwaving about what I thought would be good passband, ripple, etc. But who cares? What I really want is for |L y - F|^2 to converge quickly and for A,d,k to be "most" representative of the data.

This takes me down the road that I should be thinking of filters with L^t F = 0 so that filter doesnt change A,d,k.... It also makes me think that L could profitably help as a preconditioner. This seems like a different design optimization.

I'm also thinking it is nothing new and probably written up in textbooks everywhere. Any pointers to them would be appreciated, thanks.


texts: Causal LTI system having exponential input

texts:
I know that for an LTI system having complex exponential input,
i.e, $x(t)=\exp(j w_o t)$ & $h(t) \to $ LTI System ;
then, its output { $y(t) \} =M \exp(j w_o t + \phi)$ , where $M= |H(j w)|_{|w= w_o}$ & $\phi = \angle \{ H(j w)\} _{|w= w_o}$

But for a causal LTI system having exponential input,
i.e, $x(t)=\exp(a t)$ & $h(t) \to $ Causal LTI System ;
Is it true to say its output $\{ y(t)\}$ is given by:
 $y(t)  =K \exp(a t )$ , where $K= H(s)_{|s=a} \quad$  ?


texts: Measuring how "non-noisy" a sound signal is

texts: Given a sound clip (whose source is an FM radio signal) I'd like to measure how different from background noise (the "in-between channels" noise) the sound clip is, so that I can compare the "non-noisiness" (and thus the quality) of two radio signals. Is there a good way to do this?


texts: Deconvolving the measurement system from the system response

texts:
I'm preforming a system identification measurement, and would like to deconvolve the measurement system's response from the measured response.

The set up of the system is as follows; the signal is generated on a PC, the signal is then sent to the sound card, which in turn forwards the signal to the power amplifier, this drives a loudspeaker.  The loudspeaker propagates the sound energy which is captured by a microphone, and back into the sound card.  The recorded signal is then processed by the computer.

The signal from the sound card is also fed straight back to itself, which allows me to measure the response of the pc -> DAC -> ADC -> pc response.

I'd like to deconvolve this 'measurement chain' from the overall 'system chain'.  What methods could I use for this task?


texts: Image processing books Gonzalez et al. vs Sonka et al?

texts:
I don't know where else to post this question so here it goes. I'm trying to buy a DIP book and the best two right now are:

Image Processing, Analysis, and Machine Vision by Roger Boyle, Milan Sonka, Vaclav Hlavac 4th ed.
Digital Image Processing by Rafael C. Gonzalez, Richard E. Woods 4th ed.

They're both pricey books. So if anyone has read these books and can make a recommendation or a pros/cons comparisons of the books, I'd appreciate it.

I've already bought the "Digital Signal Processing: A Practical Guide for Engineers and Scientists by Steven W. Smith".


texts: How does Delta-Sigma Modulation convert a 1-bit signal to higher resolution signal?

texts:
I've looked at various resources about Delta-Sigma modulators, and I find them fascinating.  However, I'm confused on a particular point and and am looking for an explanation.

In my understanding, the Delta-Sigma modulators operates by sampling at a very high sampling rate (much higher than the ulitmate rate) an input analog signal.  During this process the digitized signal converted back to into an analog waveform and fed back into the input signal (possibly a number of times) which has the effect of high pass filtering the quantization noise -- so-called 'Noise Shaping'.

So far, so good.  To me, this makes sense.  However, the next stage, I find confusing.

The resulting digitized signal has been quantized using only single bit, and somehow, a digital filter is then applied to this signal which does two things:

(1) Filters out the shaped quantization noise and
(2) Re-quantizes the signal to a much higher resolution (e.g. 1-bit => 24-bit)

The resulting signal is then decimated to the desired sampling rate.

My question is about the re-quantization.  How is this performed?  Any pointers appreciated.  If you wish to poke any other holes in my understanding of DSM, that's also fine.  Cheers!


texts:
I have a continuous system described as:
$$y(t)=2x(t)+0.5x(t-2)+0.25x(t+2)$$
and I'm trying to understand why the system is stable. I know that a system is unstable when you provide an input and as a result the output is not bounded. So, when I looked at the system I set $t=\infty$ and assumed that the output goes to infinity because I'd have
$$ y(\infty)=2x(\infty)+0.5x(\infty-2)+0.25x(\infty+2)$$
making it unbounded. Looking at this equation, isn't the output going to infinity resulting in an unstable system? The solutions state it's stable and I can't see why. Thanks


texts: In the context of DFT, Where Does the Nyquist Frequency Sample Belong In a Double Sided Frequency Spectrum (Positive / Negative Side)?

texts:
If we have an even number of data points $N$, after DFT in MATLAB, the output has the order:
$$(\text{DC}, f_1, f_2, \ldots, f_{N/2-1}, f_\text{Nyq}, -f_{N/2-1}, -f_{N/2-2}, \ldots, -f_1)$$
For real signals, the first output corresponding to $k$=0, is real and so is the Nyquist frequency. After that numbers are complex conjugates.
If we are interested in a single sided spectrum, the Nyquist frequency is shown on the positive side.
However, when a double-sided frequency spectrum is plotted, many authors put the Nyquist frequency on the negative side.
Some software like OriginPro, follow the opposite. Is there a fundamentally correct way or is it just a convention i.e.,
$$
\text { If } N \text { is even, } \quad k\quad\text { takes: }-\frac{N}{2}, \ldots,-1,0,1, \ldots, \frac{N}{2}-1
$$
Alternatively,
$$
\text { If } N \text { is even, } \quad k \text { takes: } -\frac{N}{2}-1, \ldots,-1,0,1, \ldots, \frac{N}{2}$$
where $k$ is the DFT index vector, which is used to construct the frequency axis as
$$\text {Frequency axis}=k/ N\Delta t$$
where $\Delta t$ is the sampling interval.
Many people say it is just a convention and both are correct. Thanks.


texts:
Let's suppose i have a sensor that measures the output of a process and this cause a not negligible noise that affects my signal. My goal is to analyze the process signal in order to find faults. How can i deal with noise ? Is it possible to apply the fourier transform to the signal with noise or i have to filter it in advance? Maybe is possible to deduce which frequencies does the noise affect...

If i have to filter the noise, which are the techniques to do it ? Thanks


texts:
In the area of broadband spectrum sensing, there are multiple methods used to detect energy. Two common methods are:

Energy detection in the time domain (matlab example)
Spectrogram (often defined by a time and frequency resolution [Fs/NFFT], overlap percentage, and number of noncoherent averages [NAVG])

When doing energy detection in the time domain, a polyphase channelizer is often used make smaller subchannels of a large chunk of spectrum.

Questions:

How do you compare the sensitivity (ie probability of detection) between the two methods?
Is it accurate to say that the difference between the detection SNRs would be: $$ \Delta \text{SNR} = 10\log_{10}(\Delta f_1 / \Delta f_2) + 10\log_{10}\sqrt {N_{AVG}}$$ where $$ \Delta f = 1/F_s $$ for option 1 and $$ \Delta f = F_s/N_\text{FFT}  $$ for option 2?


texts: Is there a software/application that shows the FM base stations nearby us?

texts: I want to collect a list of FM base stations near me. I have found apps/software for finding LTE base station but found no success in finding an app for FM base station.


texts: binary matrices with directional non-zero elements

texts:
This question is in the context of image processing.

Assume I have a matrix of size 64 by 64. Assume even rows are all zero and odd rows are all one. I want to rotate the matrix d degrees. For example if I rotate it 90 degrees, then I will have a matrix with even columns equal to zero and odd columns equal to one.
Is there an algorithm to do that easily? Given the matrix and d I want to rotate the elements. Matlab's imrotate won't work since after rotation I get some zeros at the corners. I know one way is to use repmat, rotate it and then crop it in the center. But that does not seem straightforward. Any smart ideas? For example, is it possible to define a formula to give the values of each element given d (and the fact that the matrix was originally 1 0 1 0...?


texts: How to derive Q from r in this filter definition?

texts:
I am trying to build a resonance filter with real-time control of both center frequency and Q. I've come across what looks to me like a suitable "recipe" for a second-order IIR filter that might do the job, at http://www.music.mcgill.ca/~gary/307/week2/filters.html (toward the bottom of the page, quoted here with some loss of formatting):

$$y[n] = x[n] - a_1 y[n-1] - a_2 y[n-2]$$

where $a_1 = -2r \cos(2 \pi f_0 T) \\
      a_2 = r^2 \\
      f_0 = \text{resonant frequency} \\
      T = \text{sampling period} $

or more "refined:"

$$y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2]$$
where $ b_0 = \frac{1 - r^2}{2} \\
    b_1 = 0 \\
    b_2 = -b_0 $

The text says of the parameter $r$ that "the closer $r$ is to 1.0, the narrower the bandwidth of the resonance peak." This sounds to me like the definition of $Q$, but I suspect that it's not identical.

My question: What is the relationship between $r$ above and the more usual filter parameter $Q$?

I am just beginning to learn (self-teach) how to build filters, and I am still unclear on some of the concepts. Thanks in advance for any help.

See also this question and its answer.


texts: Maximum Likelihood carrier phase estimation of $A\cos(2\pi f_ct)$

texts:
This question refers to the example 5.2-1 of the book "Digital Communications, 5th Edition by John Proakis, Masoud Salehi" (example 6.2-1 in the 4th edition of the book) which asks for the estimate carrier phase $\phi$ that maximizes the log-likelihood function and equivalently the likelihood function.

The received signal is:

$r(t)=A\cos(2\pi f_ct+\phi)+n(t)$

where $A$ is the amplitude of the transmitted signal (known), $f_c$ is the carrier frequency (known) and $n(t) $ is the AWGN noise added.

The log-likelihood function is given as:

$\Lambda_L(\phi)=\frac{2A}{N_0}\int_{T_0}r(t)\cos(2\pi f_ct+\phi)dt$

where $\frac{N_0}{2}$ is the power spectral density of AWGN (known constant for all frequencies)
In order to get the maximum value of this we require:

$\frac{d\Lambda_L(\phi)}{d\phi}=0 \\
\Rightarrow \int_{T_0}r(t)\sin(2\pi f_ct+\hat{\phi}_{ML})dt=0 \\
\Rightarrow \hat{\phi}_{ML}=-\tan^{-1}\Bigg[\frac{\int_{T_0}r(t)\sin(2\pi f_ct)dt}{\int_{T_0}r(t)\cos(2\pi f_ct)dt}\Bigg]$

Everything OK until the last line. I don't understand how this equation is solved.


texts: How to deal with audio with different length in Audio Classification

texts:
Currently, I am trying to work with the Dataset UrbanSound8K to try some Audio classification. And I got stuck in the preprocessing step already.

Since the audios are of different lengths, like 4 seconds or 0.3 seconds, I found it impossible to directly pass into the whitening algorithms like PCA even after Feature Extraction, using mel-spectrogram/ MFCC.

So my question is what I can do under such circumstance. I was wondering about zero-padding at the end of the shorter sequence. But it seems not working and not going to yield a nice result.

I saw some people using MFCC and summarizing the MFCCs along the time-axis, like mean, variance, kurtosis, skewness.... I think that would work in this case but I just wonder if there are any other ways to do so.


texts:
I want to calculate the closed-form analytical expression for DFT coefficients of the following problem

  $$x[n]=\left\{
 \begin{array}{ll}
  1  & \mbox{if } 0\leq n \leq M-1 \\
  0 & \mbox{if }M \leq n \leq L-1
 \end{array}
\right.$$
  Write out the closed-form analityical expression for its DFT coefficients $X[k]$.

I always get the wrong solution... My attempted solution is: $$\frac{1 - e^{-j\pi k \frac{M}{L}} }{ 1 - e^{-j\pi k \frac{1}{L}} }$$

What is the correct solution and why?


texts: What are typical audio waveform synchronization methods?

texts:
I have two roughly similar (not identical) audio files, both being recordings of the same audio source material, but made independently of each other.
I need to calculate the time offset between each file, for instance one recording would be started at time A and the other at time B. By inspecting and processing the audio wave forms in each of the files I need to calculate the offset (A - B).
What are typical and effective ways of achieving this?
It can be assumed that both files are mono and that the sampling rate is consistent between the two files.


texts:
I am trying to frequency modulate a cosine signal(for assignment). I have written the following code but not getting the required output.

dt = 0.0002;
t = 0:dt:0.05;
fm = 50;
fc = 1000;
beta1 = 1;
beta2 = 2;

%making sine and m(t)
m_t = 5.*cosd(2.*pi.*fm.*t);
sin_fm = sind(2.*pi.*fm.*t); %integral of m_t. A_m adjusted in beta

%fmmod signal
fmmod1 = cosd(2.*pi.*fc.*t + beta1.*sin_fm);
fmmod2 = cosd(2.*pi.*fc.*t + beta2.*sin_fm);

%frequency_deviation
kf_1 = (beta1*fm)/5; %5 is amplitude of message signal(Am)
kf_2 = (beta2*fm)/5;

%frequency modulated signal of matlab
fmmod1_m = fmmod(m_t,fc,1/dt,kf_1);
fmmod2_m = fmmod(m_t,fc,1/dt,kf_2);

figure(1);
plot(t,fmmod1);
title('Without using fmmod beta=3');

figure(2);
plot(t,fmmod2);
title('Without using fmmod beta=5');

figure(3);
plot(t,fmmod1_m);
title('With using fmmod beta=3');

figure(4);
plot(t,fmmod2_m);
title('With using fmmod beta=5');

I have used sind in place of sin in the code because for earlier code of Amplitude Modulation sin was not giving the required waveform but when I used sind it was giving the required waveform.

AM code

dt1 = 0.0001;
fc1 = 500*1000;
fm1 = 5*1000;
t1 = 0:dt1:0.1;

%making cosine for message and carrier for higher frequency
cos_fm1 = cosd(2*pi*fm1.*t1);
cos_fc1 = cosd(2*pi*fc1.*t1);

%making Udsb_am signal for higher frequecny
Udsb_am2h = 2.*cos_fc1 + cos_fc1.*cos_fm1;
Udsb_am1h = 1.*cos_fc1 + cos_fc1.*cos_fm1;
Udsb_am_half_h = (0.5).*cos_fc1 + cos_fc1.*cos_fm1;

figure(1);
plot(t1,Udsb_am2h);
title('AM for Ac=2');
figure(2);
plot(t1,Udsb_am1h);
title('AM for Ac=1');
figure(3);
plot(t1,Udsb_am_half_h);
title('AM for Ac=0.5');

What I want to ask from the question?
When should I use sind and sin?(I know sind is for degree and sin is for radian).Asking question about sind and sin because when I saw the code on internet they have used sin not sind.

What is wrong in the code of FM?


texts: Why are LTI filters not able to output frequencies not present at the input?

texts:
I know this a very basic question but I'm just having troubles understanding why.

Thanks in advance!


texts: Can Convolution neutral network train(learn) separately (train different times)?

texts:
I am new in Convolution neutral network(CNN). My question is, is there a way to let CNN train separately?  For example, at very beginning, CNN only need to learning hardwriting 0 and 1. After the training is good enough(0 and 1 only), now, I want CNN to learn hardwriting 2. So, since 0 and 1 had learned really good, is there a way to let CNN just learning hardwriting 2 only instead of learning everything (0 1 2) from scratch ??

You know, just like human brain, after you already know how to distinguish 0 and 1,  we just need to learn 2, then we do can distinguish 0,1,2. instead of clean our memory first,then learning 0, 1,2 at same time. Because that efficiency is so low. Please advise me. Thanks.


texts: Position of poles and Stability in $z$ domain

texts: We know in Laplace Transform, if the poles lie on the left of $j\omega$ axis, we can say the system is stable. Similarly can we comment on the stability based on poles position in $\mathcal Z$-transform domain?


texts: Why do we say that spectral efficiency is measured in [bit/s/Hz], not just [bit]?

texts:
The link spectral efficiency of a digital communication system is measured in bit/s/Hz.

We have that  1 Hz = 1 s$^{-1}$.

So, $\left((1 \,\text{bit}) /1  \text{s} \right)\, / \left({\frac{1}{1 \text{s}}} \right) = 1 \, \text{bit}$.

Thus, why don't we just say then that spectral efficiency if measured in bits, but wherever I read it's bit/s/Hz? Is there some special meaning to this?


texts: Is the periodogram squared-magnitude DFT or squared-average DFT?

texts:
How is the periodogram defined? I have been searching the internet a lot, and I find two different definitions:

periodogram is the squared magnitude of the DFT
periodogram is the squared average of the DFT

which one is correct?


texts: How can I implement a phase detection algorithm with my amplitude detection?

texts:
I have been recently working on a signal detection problem in Matlab. At the minute I am detecting signals based on their amplitude, using an m-of-n detector. I am interested in improving the false alarm rate by implementing a phase detector in series with the amplitude detector.
I am very new to the field of signal processing and I'm not sure how to go about this. Here's how I think I'd go about it:
If I complex baseband the signal, that outputs a phase. I can then calculate the the autocorrelation function of the phase. If there is no autocorrelation (i.e., the ACF is close to a delta function), then it's likely that the signal is random noise. If there is some correlation for non-zero lags, then the signal correlates with itself and so must not be completely random, so must be a signal.
If this algorithm sounds like it will work, is there an easy way to combine the two methods? What do I do if the amplitude detector makes a detection but the phase doesn't, or vice versa?
As I said I'm new to this field so if this sounds like nonsense, or if there's something I'm missing, please let me know.


texts:
One of the problems of my home assignment was this.

If i use the definition to solve this(I'm not seeing other way) i have an expression like $\sum\sum nx[n]e^{-j\omega n}$ on the numerator and $\sum\sum x[n]e^{-j\omega n}$ on denominator. How do i simplify those expressions? The double sum is making me confuse.

Thanks.


texts: Given the Power Spectral Density of a random process determin if it is stationary in a wide sense

texts:
I am studying a Communication Systems course and have trouble figuring out how to tell and prove if a random process is stationary in a wide sense while only having it's power spectral density.
some examples of the problems I am tackling are,
S(f) = sin2πf/(2πf)
S(f) = cos(2πf)^2
S(f) = (cos2πf)^2
The only thing I understand proves a Random Process is stationary in a wide sense is that it's autocorralation function is subject only to the time shift τ and not time, but I haven't seen any rules on how it can be determined from the Power Spectral Density.


texts:
I am trying to design a high pass filter to remove baseline drift from an ECG signal.
the baseline drift is of very low frequency like 0.3Hz or so with an amplitude of 25% of the ECG signal. But the bandwidth of ECG signal itself is 0.5Hz to 150Hz. I have been trying fdatoolbox in matlab to design the HPF but I was not able to remove the baseline drift accurately.Please help on how can i achieve a low order filter to eliminate baseline drift.
I have got it using higher order filters of order above 1000. But I'm looking for a way if I can get lower order filter,like order of less than 30, to do the same.Thanks for your help..


texts: Which order to perform downsampling and filtering?

texts:
I have EEG data recorded with 128Hz sampling rate.
As my goal is to reduce the amount of data (and maybe noise), I want to downsample the data to 64Hz (I am only interested in range 0.5 - 30Hz).

I would perform downsampling (64Hz) and use a bandpassfilter (0.5 - 30Hz).

Now my question is, what is the best order to perform those two steps and why?
Which impacts will these methods have to the quality of my data?


texts:
I am trying to compare several PSD methods to determine the most suitable method for my application. Does anyone know a good approach to start with? I am relatively new to this field of topic.
I could possibly compare the results of the psd methods, the computational complexity and the resolution of the results.


texts: How to design a *desired signal* for Lyapunov stability analysis?

texts:
How to design a desired signal for Lyapunov stability analysis?

Like given in this paper:

http://journals.tubitak.gov.tr/elektrik/issues/elk-15-23-3/elk-23-3-8-1212-29.pdf


texts: How to extract the Freezing Index from Accelerometer data?

texts:
I'm working on a project to detect Freezing of Gait episodes for Parkinson's Disease patients, and according to this paper and other ones, extracting the freezing index which is "The power in the “freeze” band (3-
8Hz) divided by the power in the locomotor band (0.5-3Hz)" will be a good feature.

Here is the explanation: One standard feature which is extracted from the raw signals is the Freezing Index (FI), defined as the ratio between the power contained in the so-called freezing and locomotion frequency bands (3-8 Hz and 0.5-3 Hz respectively). This feature is convenient since it requires only FFT-computation.

But I don't know how to implement it in Python.

I have a dataframe like this:

Then, I did something like this to extract features from the sensors time-series data:

win_size=200
step_size=40
for col in df.columns:
    if col != 'Label':
        df[col + '_avg'] = df[col].rolling(win_size).mean()[step_size - 1::step_size]

Now, I would like to extract the Freezing Index, How can I do this? and I need someone to explain it to me because I don't fully understand it!


texts: Analysing DAC Spectra: Transient Noise Analysis

texts:
I am working with a new Digital-to-Analog Converter (DAC) design in simulation and I'm trying to analyse the output. The device takes in an ideal 14-bit digital representation of a sine wave and outputs through an ideal Butterworth filter (a̶n̶t̶i̶-̶a̶l̶i̶a̶s̶i̶n̶g̶  anti-imaging).
In the simple analyses, I have been setting my input sinusoidal frequency such that it is coherent, and my sample capture length set up to collect 8196 points. I'm then running this through a simple MATLAB script that windows the simulation data and calculates the Welch periodogram using pwelch(), to reduce spectral leakage and variance respectively. Once I've done that, I'm measuring the characteristics of the device such as the SNR, SINAD (SNDR), THD and so on.
The next step in my analysis is slightly more complicated. I have run 5 transient noise simulations which result in 5 sets of output data, the only differences being the noise seen on the transient signal. How can I combine these spectra to get a more accurate picture of where the noise floor is in the device?
Assuming the noise is random, can I approach the situation similar to Bartlett's method but averaging across all 5 FFTs and creating a periodogram of the result?


texts: An Impulse in Complex Envelope (i.e. baseband) Notation?

texts:
In complex envelope notation, a modulated RF signal, say v(t), is represented in complex envelope notation as v(t) = Re{ g(t)e^(jwct)}, where wc is the carrier angular freq, g(t) = x(t) + jy(t), and g(t) is the complex baseband signal.  Since x(t) and y(t) are inputs to an in-phase (cos) and quadrature phase (sine) mixers, how would I represent v(t) in baseband g(t) notation, where v(t) = an impulse?  Since a mixer is technically nonlinear, and therefore NOT an LTI system, isn't this impossible to do in a mathematically correct way?
The issue is we have a MATLAB RF simulator for PHY signals, whereas the simulator is designed to work only with RF signals represented in baseband notation (and fc is metadata), and I want to send impulses through the RF channel to clearly see the fast fading reflections.  My current method is to use impulses for x(t), with y(t) = 0, but not sure if this translates correctly, or if this even makes sense.


texts: FIR filter design by the Fourier transform method

texts:
I am having some problems understanding how the Fourier transform method is used to determine the FIR filter.

As far as i have understood, you start by using the ideal impulse reponse for the specific filter design you want to have.  Those are given in this table

As the the ideal impulse reponse contains infinite impulses, we truncate the infinite to finite a number by $m$, (which is confusing how this is done).

But my question is how do predetermine the order of the filter, without having to try randomly, using this method?


texts: Trellis-converging offset and the min-sum algorithm

texts:
I am going through a paper on the design of interleavers for Turbo codes and I am a little bit confused about how to calculate what they call Trellis-Converging offset. They define such parameter of a convolutional code as

  Trellis-converging offset is defined as a doscrete random variable ($W$) on positive integers. The outcome value $w$ corresponds to the event that the trellis converges in $w$ number of symbols in a convolutional code decoder using the min-sum algortihm.

The definition is pretty clear by itself, but I am confused with what they mean by convergence. As far as I know, the min-sum algorithm operates using a fixed set of hidden states on the factor graph passing messages between them, and so if no fixed number is set, the algorithm would need to get the full stream of bits to decode the sequence. However, they say that such algorithm converges for a number of symbols that is a random variable. I have been exhaustively searching in the internet and I have not found information about this TCO, so I was wondering if someone could clarify me this converging criteria or give some reference about the topic.


texts:
I am trying to extract the decay properties of an oscillating time history. The obvious method for doing this is to use a Hilbert transform. The figure shows an example of a theoretical decay and the envelope calculated using the Hilbert transform.

My problem is that the envelope has oscillations (Gibbs oscillations) due to the fact that there is a discontinuity in the sudden start of the data.

The extracted decay rate, below, which should be a constant in this theoretical case, is unacceptably corrupted by the oscillations.

I have considered filtering the envelope with a low pass filter but as the oscillations are the same order as the original frequency I will spoil the decay data. Alternatively I could create a mirror image of the data and thus remove the starting step. However, this requires finding the starting phase of the time history. Could these methods work or can you suggest another method?

Edit
As suggested by Olli Niemitalo I have had another look at filtering the envelope signal. The filtered envelope is added here

The filtered envelope signal now has no oscillations but has dropped down from the peaks. This is fine for my application were I need the rate of decay. The next figure compares the decay calculated from the filtered envelope to the unfiltered envelope

This is success! Thanks to Olli.

Thanks


texts:
I am trying to include the ability to have sub-clock delays on an FPGA system's data bus without the need of fractional delay filters.  I am using a Xilinx V6 part, and could utilize its IODELAYE1 instance, but the max delay it can create is less than half a period, and the parts cannot be cascaded.

Is there theory for how to to do sub clock delays of data buses within a system?


texts: Unknown symbol/expression in text about adaptive filters (cst)

texts:
I am currently reading a chapter about adaptive filters from the Springer Handbook of Speech Processing.

In a formulation of the variable stepsize normalized least mean squares (VSS-NLMS)-algorithm, i found an expression that I do not know yet, and searching the internet did so far not yield results i can make sense of.

The algorithm is expressed as

$$\mathbf{\hat{h}}(k) = \mathbf{\hat{h}}(k-1) + \mu_\text{NPVSS}(k)\mathbf{s}(k)e(k),$$

where $\mathbf{\hat{h}}(k)$ is the current estimate of the desired filter impulse response, $\mathbf{s}(k)$ the source signal and $e(k)=x(k)-\hat{x}(k)$ the error difference between filter output and actual output (using the real $\mathbf{h}(k)$).
The stepsize for this algorithm is defined as
$$\mu_\text{NPVSS}(k) = \frac{1}{\mathbf{s}^\text{T}(k)\mathbf{s}(k)} \left[ 1-\frac{\sigma_b}{\sigma_e(k)} \right],$$

with $\sigma_b$ and $\sigma_e$ denoting the square root of the variance/power of a noise process $b(k)$ uncorrelated to the source signal and of the error difference, respectively.

For a practical implementation it is suggested to compute the fraction $[\mathbf{s}^\text{T}(k)\mathbf{s}(k)]^{-1}$ as $[\delta + \mathbf{s}^\text{T}(k)\mathbf{s}(k)]^{-1}$ instead.

In this context, the regularization parameter is defined as
$$\delta = \text{cst} \cdot \sigma_s^2.$$

I am not too familiar with the concept of regularization, so I have no clue what the term $\text{cst}$ means in this algorithmic formulation.

What does it mean? Any pointers or literature recommendations would be very much appreciated.


texts:
I have a system where I sample some data periodically (every 10usec). The shape of the data is triangle, in other words it linearly increases and decreases in time. (Both theory and practice are in agreement, I plot the ADC output and I clearly see linear increase and decrease although there is some non-linearity but let's ignore that for the sake of discussion) However the frequency of the signal is not fixed. It may increase very fast or very slow. (narrow or wide triangle)

I am trying to locate the true peak point of the data with 500nsec accuracy, however my ADC samples every 10usec. One item on my favor is that the data values increase and decrease in sequence. I get worst case 50 samples, best case 1000 samples. (in worst case, I will have 25 samples in increasing order and 25 in decreasing order)

How can I approach this issue using DSP techniques? One idea I had was to calculate peak point from rise and fall angles. (I basically extend the left and right side of the triangle and look at where those lines cross each other). I cannot think of any other way. I welcome suggestions.


texts: Couting exact number of mutiplications in Separable filters with Zero Padding

texts:
I'm convolving an Image $I$ of dimension $2\times 3$ with a kernel of size $2 \times 2$ and counting the exact number of multiplications.My image and kernel are given as follows
$$
I=\left[\begin{array}{ccc}
2 & 0 & 1 \\
1 & -1 & 2
\end{array}\right] F=\left[\begin{array}{cc}
1 & -1 \\
1 & -1
\end{array}\right]
$$
Now, as there is zero padding according to Tensorflow's convention I am padding with zeros on bottom and right. Now My image and filter after flipping becomes,
\begin{align*}
\mathbb{F'} & = \begin{bmatrix}
-1 & 1 \\
-1 & 1
\end{bmatrix}
\mathbb{I'} = \begin{bmatrix}
2 & 0 & 1 & 0 \\
1 & -1 & 2 & 0 \\
0 & 0 & 0 & 0
\end{bmatrix}
\end{align*}
Now doing correlation b/w $\mathbb{F}$ and $\mathbb{I}$ we get
\begin{align*}
  (F' * I') =  \begin{bmatrix}
-1*2+1*0+-1*1+1*-1 & -1*0+1*1+-1*-1+1-2 & -1*1+1*0+-1*2+1*0 \\
-1*1+1*-1+-1*0+1*0 & -1*-1+1*2+-1*0+1*0 & 2*-1+1*0+-1*0+1*0
\end{bmatrix}
\end{align*}
\begin{align*}
    (F' * I') &= \begin{bmatrix}
-4 & 4 & -3\\
-2 & 3 & -2
\end{bmatrix}
\end{align*}
The total number of multiplications here are, 4 * 3 * 2 = 24.
We know that if a kernel is separable the $rank$ of kernel will be 1 . Here as the second column of the filter $\mathbb{F}$ is a scalar multiple of the first column, the $rank(\mathbb{F})=1$.
\begin{align*}
    \mathbb{F} &= F_1\cdot F_2 \\
    &=\begin{bmatrix}
    1\\
    1
    \end{bmatrix} \times
    \begin{bmatrix}
    1 & -1
    \end{bmatrix}
\end{align*}
\begin{align*}
    \mathbb{I'} &= \begin{bmatrix}
2 & 0 & 1 & 0 \\
1 & -1 & 2 & 0 \\
0 & 0 & 0 & 0
\end{bmatrix},
    F_1 = \begin{bmatrix}
    1 \\
    1
    \end{bmatrix},
    F_2=\begin{bmatrix}
    1 & -1
    \end{bmatrix}
    \\
    (F_1* I')&=\left[\begin{array}{cccc}
3 & -1 & 3 & 0\\
1 & -1 & 2 & 0
\end{array}\right] \\
(F_2*(F_1*I')) &= \left[\begin{array}{ccc}
-4 & 4 & -3\\
-2 & 3 & -2
\end{array}\right]
\end{align*}
Here the total number of multiplications are (2 * 4 * 2) + (322) = 28
But if the kernel is separable the total number of multiplications should be less than the total number of mutliplications if we do direct convolution. Where am I going wrong?
Direct $2 \mathrm{D}$ convolution is $O\left(M_{1} N_{1} M_{2} N_{2}\right)$
while two $1 \mathrm{D}$ convolutions is $O\left(M_{1} N_{1}\left(M_{2}+N_{2}\right)\right)$


texts: Feedback command for MIMO systems in matlab

texts:
I have a state space model with 12 outputs and 4 inputs. I want to close the loop of 8th output with 3 input. then in outer loop i want to close outer loop for 1 output and 3 input again (Quadrotor Pitch and X position states).
I have close first loop using this command
S1=STATE_SPACE_MODEL(8,3);
S_CL=feedback(S1*PID,1);

But, now i m confused how do it close the outer loop. Please help me out in this!!
this has to be done using feedback command.


texts: How to evaluate performance of a simulated system

texts:
I have created a database (DB) of amplitde vs time delay of UWB pulses originating from various (known) points in a certain (known) "region" in a space. This DB also stores region name, since there are several regions . This is "training phase".
the folowup part is called "localization phase" -- I assume that there is a signal coming from one of the regions ..I don't know which, I changed parameters of simulation and now know this new comer's amplitude vs delay profile. I will cross-correlate this new signal with the ones stored in my DB in #1 above. If any of the signals matches, I can then say that the signal was from region 1 or 2 or 3 .. etc.
My problem: How do I evaluate this system? That it fared better or worse? If I had resources, I could probably buy or borrow a network analyzer, antenna, etc, and measure the delay profile and that would give me something to compare against.

Now without that option, what are my (other) options?

Gurus, please advise.
Thank you.


texts: Is the Hilbert transform computable "losslessly"?

texts:
The Hilbert transform is a procedure that can be used to perform a 90° phase shift on a signal.

Is the computation of (a discretized) Hilbert transform a "lossless" process? I.e. it merely phase shifts the signal, but does not alter it in any other way.

The reason for asking is that the following picture here looks like it introduces slight variation in the Hilbert transformed signal.

http://www.originlab.com/doc/Origin-Help/Hilbert-Transform


texts: How to calculate SNR before and after filtering for noise cancellation in an ecg signal in matlab?

texts:
I am working on noise cancellation in ecg signal and I want to calculate the Signal to noise ratio(SNR) before applying filter to the signal and after applying filter i.e., snr_before and snr_after. I used the following commands:

snr_before = mean(input_or_desired_signal.^2)/mean(noise_signal.^2);
residual_noise = desired_signal - filtered_output_signal;
snr_after = mean( desired_signal .^ 2 ) / mean( residual_noise .^ 2 );

The code works as no error is generated.
But only the box is displayed when I plot d same in each case and no graph is generated..
Pls help
Thanks in advance..


texts: Direct Inverse Filtering an image, spatially convolved with a blurring function without any noise, doees not give my original image back?

texts:
I am trying to learn about inverse filtering of images and am trying it out on a test image in MATLAB.

The test image is an 8x8 checkerboard, with 64x64 pixels, generated by the MATLAB code:

g = checkerboard(8);

I am using a 3x3 averaging filter, h, and convolving it with the test image.

I am not adding noise and I want to see the results of direct inverse filtering the distorted image.

I am applying the deconvolution function in the frequency domain by:

1) Padding the 3x3 averaging mask h to generate a 64x64 mask hp.

2) Centering the coefficients (0.111...) of hp to maintain symmetry.

3) Multiplying hp by (-1)^x+y (since I want the fft of the mask to give me a centered Fourier transform, hpc.

4) Taking the FFt of the padded, centered mask and multiplying (-1)^(u+v) (since I centered the mask in the spatial domain in step 2) to give me H.

5) Taking the distorted, and frequency shifted test image and dividing it by H, (ie. Y = G./H). (Note: I am not padding the image since the mask was padded, and wraparound is not a huge concern here)

6) Obtaining the inverse filtered image inv = (-1)^x+y * real(ifft2(Y))

My question: Why can't I obtain the original (or something very close to it)?


texts: Uncertainty principle - Duration bandwidth principle

texts:
In one of NPTEL courses about time-frequency analysis, the professor said that the duration bandwidth principle is $\sigma_t^2 \sigma_\omega^2 \ge \frac{1}{4}$.

He added that the formula making use of time resolution and frequency resolution is a false one. The time and frequency resolutions correspond  to the time "distance" between two samples of the signal and frequency resolution corresponds to the "distance" between two samples in the frequency domain (two successive samples of the Fourier transform of the time domain signal). Can anyone, please, clarify this to me as I get to see the one using time and frequency resolutions in so many papers?

NB: Here is the link to the course I'm talking about http://nptel.ac.in/courses/103106114/


texts:
I have 2 images of the same scene, captured by different cameras. So they a bit displaced. As far as I can understood, phases of FFT contains info about image details. My goal is to take intensity from first image, details from second and produce enhanced highly detailed image with proper intensity.

I tried to solve problem as is - apply FFT for both images, take magnitudes from first image and phases from second and apply iFFT. But I didn't get what I want. Image looked corrupted, see below.

Then I tried to compute Phase correlation , align images and
apply the algorithm described above again. Result was the same.

Is it possible to do what I want? May be I need to do some calculations in frequency domain, like phase shifting? Matlab code is highly appreciated!


texts: Design Filter with Arbitrary Magnitude and Phase Response in Matlab

texts: Does anyone know of a way to design a filter with arbitrary magnitude and phase response in Matlab?


texts: What is the significance of first zero upcrossing of a cross correlation?

texts:
What is the significance of the first zero up-crossing of a cross-correlation? The signals being correlated are cosines signals with some phase difference.

EDIT 1 : Cross-correlation plot attached. The signals are cosines with phase difference and may be with some low frequency components.

Sampling frequency of the measured signal is taken as 50Hz

Cross correlation obtained using:

Cxx=fftshift(ifft(fft(x,N).*conj(fft(y,N))))/(norm(x) * norm(y));


texts: Implementation of Hilbert Spectrum in MATLAB

texts: I am trying to implement this equation  on Matlab,  where i have wj values in the form of an array, cumtrapz command of MATLAB implements this, but its confusing which way is exact. such as cumtrapz(t,wj) or simply cumtrapz(wj) or manually doing it? Because results of all these though related, are different.


texts: On the de-noising of time-frequency images

texts:
I am wondering what techniques might be available for 'de-noising' the following example time-frequency image that was created using Welch's method. The following plot was created from a robotic sensor. (This is NOT a color image - it is a greyscale image - colors added for visual purposes only).

Goal:

My goal ultimately is to estimate the pulse-spacings that you see here, should such pulses exist. This might be somewhat of a chicken and egg, so to this end, I ask myself, "Do pulses of this rep-rate +/- 10% exist?", and go about detecting them. What you are seeing here is the signal (pulses), but along with other un-wanted interference. However as Emre has suggested, they have structure, albeit in the Time-Frequency space. Do time-frequency filters as such exist?

I would strongly like to see image processing solutions applied here, but am open to any solution.

Thus: The goal is to remove all the high intensity signals except the repetitive pulses (found near index 300 on the y-axis) as can be seen. All the other high intensity signals can be regarded as 'interference'.

Assumptions you may make:

You may assume that you roughly know the pulse lengths that you are seeing here. (Let us say, within +/- 10%). Put another way, you have decided to look for pulses of this length. (+/-)
You may assume that you also roughly know the rep rates of the pulses, (again, let us say +/- 10%).
Unfortunately you do not know their frequency any more accurately. That is to say, in this image the pulses are at 300, but they could have just as easily been at 100, or 50, or 489, or whatever. However, the good news is that those frequencies shown here are very close to one another, on the order of say, 10's of Hz).

Some thoughts of mine:

Image processing POV:

Morphological operations have occurred to me, however I am not too familiar with those to know if they might work or not. I suppose the idea might be to 'close' and hence remove the 'bigger' stains?
Row-wize DFT operations might indicate which rows to null out, based on the rows of interest having the highest repetitive pattern, however it might not be a viable solution if the pulses are few and far between, or if image is more noisy.
Just by looking at the image, you almost want to 'reward' isolation, and 'punish' connectivity. Is there an image processing method(s) that accomplish this sort of operation? (Morphological in nature again).

What methods might help here?

Signal processing POV:

The frequency range shown here is already extremely tight so I am not sure notch filtering operations will help. Moreover, the exact frequency of the pulses shown within this tight range is not known a-priori.
By making educated guesses on the pulses of interest here, (their lengths, and repetition times) might I be able to compute the 2-dimensional DFT of my 'template', and utilize this as a 2-D cepstral-temporal filter to which I simply multiply the Welch image shown above by, and then perform an inverse 2-D DFT?
OTOH perhaps would Gabor filters be a good match here? After all, they are orientation sensitive filters, similar to our own built in V1 visual processors. How might they be exploited here?

What methods might help in this domain?

Thanks in advance.


texts: Determine the time length of audio training samples?

texts:
Currently I am trying to use MFCC feature for audio training in order to do human voice classification (can later distinguish the person who speaks).

I have problems to determine the proper time length of audio training samples? To be clear, say I have computed 12d-MFCC features over 20ms window (10ms overlapping), if the length of the audio training samples is 1s, this will result in a 12 x 100 feature vector (as there are 100 frames during the 1s time given 10ms shifting time step). In practice, what's the proper time length of audio training samples? Of course, long time length will results in long MFCC feature vector.


texts: OpenCV Haar Classifier - best type of positive images

texts: I'm currently collating images to try and develop a HAAR classifier using the OpenCV traincascade functionality. The objects I want to detect are small and numerous (bees), and I'm curious as to whether it would be better to train the image with lots of small cutouts of the bees as positive images, or to use bigger positive images with the bees present along with a chunk of the background.


texts: How can I perform sound classification based on spectrograms and waveforms?

texts:
I noticed that almost all machine learning enabled sound classification approaches use spectrograms or MFCCs as input for feature extraction.
I am new to the area and would like to know if we can also combine waveforms and spectrograms/MFCCs of sounds to classify them more accurately.
Then, is there a deep learning approach that does not require any pre-processing and relies simply on raw audio data?


texts: Finding Frequency Differences In Differing Amplitude Signals

texts:
I have two audio signals, both with supposedly matching frequencies at each time window.  What I would like to do is to verify that the frequencies are identical.  However, one was created at a higher power than the other, so simply checking for a match at each sample won't work.  I've already cross-correlated the signals, so there's no time lag between them.

Is there any way to compare the "shapes" of two waveforms of differing amplitudes?  If there's some way to drop the amplitude domain from the signal, that would work - I don't know nearly enough about DSP to know if it's easy or not.


texts: Optimal sampling rate for the HMM forward-algorithm

texts:
I have a system with a binary-state. The system state is estimated by an HMM forward-algorithm. Also, the system allows a varying sampling rate.
Considering that the system state transition takes a known time $T$, in some cases, the sampling can be reduced without a reduction in performance. How to optimize the detection of transition between the latent states and the sampling rate?


texts: Why does FFT not have an effect on my smoothed signal?

texts:
I'm playing with FFT at the moment and I try to get periods from noisy signals by recreating this example. While experimenting, I've noticed that after smoothing a quite noisy signal, the result of fft() is actually the same signal again - which is what I don't understand.

Here is a full example which can be run in an IPython Notebook (You can create a notebook here and run the code if you want).

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

figsize = (16,8)
n = 500
ls = np.linspace(0,2*np.pi, n)

x_target = np.sin(12*ls) + np.sin(52*ls)

x = np.sin(12*ls) + np.sin(52*ls) + np.random.rand(n) * 3.5
x = x - np.mean(x)

x_smooth = pd.rolling_mean(pd.DataFrame(x), 14).replace(np.nan, 0.0).as_matrix()
x_smooth = x_smooth - np.mean(x_smooth)

x_smooth = np.roll(x_smooth, -7)

# Getting shwifty and showing what we've got
plt.figure(figsize=(16,8))
plt.scatter(ls, x, s=3, c=[1.0,0.0,0.0,1.0])
plt.plot(ls, x_target, color=[1.0,0.0,0.0, 0.3])
plt.plot(ls, x_smooth)

plt.legend(["Target", "Smooth", "Noisy Data"])

# Target
x_fft = np.abs(np.fft.fft(x_target))
pd.DataFrame(x_fft).plot(figsize=figsize)

# Looks like it should
x_fft = np.abs(np.fft.fft(x))
pd.DataFrame(x_fft).plot(figsize=figsize)

# Plots the same signal?
x_fft = np.abs(np.fft.fft(x_smooth))
pd.DataFrame(x_fft).plot(figsize=figsize)

Below you find the resulting plots of this script.

Noisy data with smoothed signal and target function:

FFT of the target

FFT of the noisy data

FFT of the smoothed data

I don't really get why this is the case here. Can somebody explain this to me or am I doing something wrong here?