texts: What are the differences between these two Gabor filter functions

texts:
I need to enhance the visibility of veins on dorsal hand vein images in my project. I use two different even-symmetric Gabor filters bank improve vein visibility.

First bank consists of these gabor functions:
$$G^\mathit{e}_\mathit{mk}(x,y)=\dfrac{\gamma}{2\pi\sigma^2}\exp\Bigg\{-\frac{1}{2}\left(\dfrac{x_\mathit{\theta}+\gamma^2y_\mathit{\theta}^2}{\sigma^2}\right)\Bigg\}\times \left(\cos(2\pi f_\mathit{0}x_\mathit{\theta})-\exp(-\dfrac{\upsilon^2}{2})\right)$$

Second bank consists of these ones:

$$G^\mathit{e}_\mathit{mk}(x,y)=\exp\Bigg\{-\frac{1}{2}\left(\dfrac{x_\mathit{\theta}+\gamma^2y_\mathit{\theta}^2}{\sigma^2}\right)\Bigg\}\times \cos(2\pi f_\mathit{0}x_\mathit{\theta})$$

where $m$ is the scale index, $k$ is the orientation index, $f_\theta$ is the filter center frequency, $\sigma$ is the standard deviation (often called scale), $\gamma$ is the aspect ratio of the elliptical Gaussian envelope, $\upsilon$ is the factor determining DC response, $x_\theta=(x\cos\theta+y\sin\theta)$ and $y_\theta=(-x\sin\theta+y\cos\theta)$ are rotated versions of the $x$ and $y$ coordinates.

I have coded these filters in MATLAB, I do not have any problem on coding. But I cannot understand the underlying difference between these two gabor functions.


texts: why does butterworth IIR filter for a square pulse show ripple on edges in time domain but not the same for sine wave

texts:
I am using butterworth filter for both square pulse and sine wave. In case of square pulse, butterworth produces some sharp ripple of edges of square pulse in time domain for varying butterworth order (wn > 2). But in case of sine wave, I don't see any ripple. Could anyone please explain the reason behind that? Matlab code is following along with plots.
close all

clear all
T = 1/10;
t = linspace(0,T,1001);

s = square (2*pi*100*t);

fs=10/200; % sampling frequency
fc=0.2; % fc normalized
offset=0;
amp=3.3;
duty=50;
t=0:0.01:100;%100 seconds
sq_wav=offset+amp*square(2*pi*fs.*t,duty);
sin_wav = offset + amp*sin(2*pi*fs.*t);
disp(length(sin_wav));

N = length(sq_wav);
disp(N);
#disp(t(1:10));
t1 = [0:N-1];
disp(t1(1:10));
disp(fs);
wn = 15; % filter order

[b,a] = butter(wn, fc);

[h,w] = freqz(b,a);

out_sq = filter(b,a,sq_wav);
out_sin = filter(b,a,sin_wav);

subplot (4, 1, 1)
plot(t(1:end),sq_wav(1:end))
title("square pulse");

subplot (4, 1, 2)
plot(t(1:end),sin_wav(1:end))
title("sine pulse");

subplot (4, 1, 3)
plot(t(1:end),out_sq(1:end))
title("square pluse after butterworth");

subplot (4, 1, 4)
plot(t(1:end),out_sin(1:end))
title("sine pluse after butterworth");

Here is the output of this code:


texts: Is frequency resolution ("bin size" inverse) always directly proportional to sampling duration (window length)?

texts:
I am new to DSP and signal processing in general (digital or analogue). In implementing an FFT/Goertzel I have noticed that the resolution ("bin size" inverse) is directly proportional to the sampling duration. In other words; the tighter you want the resolution, the longer you have to sample for (regardless of sample period / how quickly or slowly you sample (above the Nyquist rate obviously)). Is there any way around this (using some technique other than FFT/Goertzel) or is this an inherent scientific principle? If so then does the principle apply to analogue filters as well?
Would a Goertzel algorithm be the best solution for detecting the presence of a particular frequency, with the best resolution (tightest filter)?


texts: problem finding 16 point DFT using two 8 point FFT (Divide and combine algorithm) MATLAB

texts:
Write MATLAB code that determines and plot the N-point Discrete Fourier Transform of x[n] defined by the following equations:

x[n]=0.5*pi*n      n=0:16

Compute and plot 16-point DFT using two 8-point FFTs and combining techniques?


texts: Why does each node in the wavelet scattering transform split into multiple paths?

texts:
Why does each node in the wavelet scattering transform split into multiple paths as in this figure from https://arxiv.org/abs/1304.6763?

I understand roughly what’s happening along a single path, but I don’t understand why each node splits into three.
Thanks!


texts:
I need to detect a known number of simple markers (chessboard patterns) in a color image efficiently.

I implemented normalized cross correlation (NCC) as a pixel shader via WebGL since I need this to work in a client application in a browser. Detection works fine and is sufficiently robust, but I need better performance. It currently takes a number of seconds to calculate the NCC image, and I'm already using 'obvious' optimizations such as cropping the image (I have a rough indication where the markers should be) and reducing the resolution. I also optimized the code itself (e.g. using loop unrolling, but that is another matter and probably belongs on SO).

The image has a resolution of 900x520 pixels and the markers are 22 square pixels. The code runs on tablet computers with Intel Atom processors.

My question is: Are there algorithms that are better at utilizing what I know about the image and the marker geometry, ideally ones that are suited for almost branchless pixel shaders? For instance, I was considering a band pass filter (a fourier transformation that only leaves the frequency of the marker, so to speak) or FNCC (JP Lewis's paper), or something that utilizes the fact that the markers are black and white.

Unfortunately, it seems these approaches are comparatively expensive since the FFT of the image isn't cheap and my markers are small compared to the image resolution. Testing a ton of algorithms blindly isn't an option, unless there's a library I'm not aware of that helps me implement this in JavaScript or as glsl pixel shaders. Even then, an indication of a more efficient algorithm would be helpful.

EDIT:
Added an image of the marker. The background is normally transparent, of course, i.e. the NCC will ignore those pixels


texts:
I am very new to this things. Sorry for probably stupid question. I don't understand what units and meaning have the values on Y-axis of Fourier Transform graph? On X-axis it is Frequency (Hz). Pretty clear. Can't find the meaning of Y.

Here is the code in MATLAB I used to produce the graph, and the image of the output:

A = 5;
t = 0:1/1000:1;
f1 = 7;
s = A * sin(2 * pi * fi * t);
F = fft(s);
figure;
plot(abs(F)), (s);


texts: Inverse $\mathcal Z$-transform of rational functions

texts:
What will be inverse $\mathcal Z$-transform for this function:

$$H(z) = \frac{\left(1+\beta z^{-1}\right)\left(1+\beta z\right)}{\left(1+\alpha z^{-1}\right)\left(1+\alpha z\right)}$$


texts:
I'm trying to implement this filter in title by following the book "Designing Audio Effect Plug-Ins in C++" By Will Pirkle

Problem I'm facing is the magnitude response error when comparing against analog model of the LP filter. Here is a plot showing the difference.

I get these are values for variables (fc=1000, fs=44100):

>> massberg_test

    w0 =  0.14248
    g1 =  0.045305
    gm =  0.70711
    wm =  6283.2
    Om =  0.071359
    Os =  0.10081
    g0 =  1.1008
    a0 =  0.14612
    a1 =  0.055508
    b1 = -0.89919
Real:
    a = [ 0.132736  0.050424 ]
    b = [ 1.00000  -0.81684 ]

I don't have access to the original paper by M. Massberg so, I can't be sure if the equations in my source book is correct (I've read that the "ac" in g1 formula is just a typo and book errata does not mention other errors for this fileter).

Q:
1. Are the equations found in Pirkle book correct?
2. If 1 is Yes then where to look the issue from?


texts:
Let there be a discrete signal that is the sum of sinusoidals and can be described by
$$s(n)=\sum_i A_i \mathrm{sin}(2\pi f_i \frac{n}{f_s}+\phi_i)$$
where $A_i, f_i, \phi_i$ are unknown but fixed parameters of my signal, $f_s$ is my sampling rate, and $n$ is the sample measured at time $t=\frac{n}{f_s}$.
Of this signal I measure the first $N$ samples ($n=1,2,3...N$) and I also measure another $N$ samples at points $n=N+20000,....2N+20000$.
I am now interested in estimating/predicting/extrapolating the value of a single sample that is exactly in the middle between my measured chunks: $s(n=N+10000)$. I don't care about all the other samples in the gap.
What is the best and most efficient way to estimate my desired sample?
Note that $N$ would be ideally between 100 to 1000 and my sampling rate $f_s=100\mathrm{MHz}$.


texts: Continuous Wavelet Transform time vector in python

texts:
I have a signal sampled at 128 Hz. I used to extract features with the spectrogram function and I decided to upgrade my algorithm and I'm trying to analyze it using  Continuous Wavelet Transform  (pywt.cwt) in python.
this function has only 2 outputs: coefficient and frequency, while spectrogram returns the time vector as well.
Similar to what is said in theory, the resulting frequency vector has varying differences between each two samples (high density at low frequencies, and vice versa). According to the uncertainty principle, the time vector also has to change its density inversely to the frequency vector, is not it? My understanding is based on the famous graph in the picture attached.

If I'm not mistaken, the time vector (that is not returned) is just the original time vector of the signal, meaning the difference between each two adjacent samples is constant.
Please help me understand this mismatch.


texts: Implementation of FFT algorithm in C Programming

texts:
Does anyone know any FFT algorithm that can be implemented in C?
I am trying to convert a built-in matlab FFT api to a function in C.

Appreciate any help!


texts:
I have a filter where the input is:
$$ x(t) = s(t) + w(t) $$
which is signal plus white Gaussian noise and the output is $y(t)$. If the impulse response of the filter is $h(t)$, then:
$$y(t) = s(t)*h(t) + w(t)*h(t)$$
I have given the following equation to calculate the SNR:
$$SNR = \frac{\int_{-\infty}^\infty|s(t)*h(t)|^2 dt}{\int_{-\infty}^\infty|w(t)*h(t)|^2 dt}$$
And in frequency domain, following the Parseval's theorem:
$$SNR = \frac{\int_{-\infty}^\infty|S(f)H(f)|^2df}{\int_{-\infty}^\infty|W(f)H(t)|^2 df}$$
I wonder how this represents the signal power and noise power at numerator and denominator, respectively? Cause it is more like calculating energy than power.


texts: How to derive the answer to this convolution problem?

texts:
I came across this below question, (which was a homework assignment question for Signal Processing class, which my friend mailed me for help solving), mulled over it for an hour and had no idea how to proceed with solving it.

Let $C(x) = A(x)B(x)$ where:

$$A(x)=\sum_{n=0}^{N_1}a(n)x^n$$
$$B(x)=\sum_{n=0}^{N_2}b(n)x^{2n}$$
$$C(x)=\sum_{n=0}^{N_3}c(n)x^n$$

Find expressions for $N_3$ and $c(n)$ as functions of $N_1$,$N_2$, the $a(n)$ and $b(n)$.

Apparently, it has something to do with convolution, as in $c(n)$ is the convolution result of $a(n)$ & $b(n)$ or something like that.
But I still can't figure out how that is. Can anybody please explain me the answer of $N_3$ and $c(n)$


texts: Fast Fourier Transform and Clustering of Time Series

texts:
I'm making a project connected with identifying dynamic of sales. That's how the piece of my database looks like:

There are free columns:

Product - present the group of product
Week - time since launch the product (week), first 26 weeks
Sales_gain - how the sales of product change by weeks

In the database there is 3302 observations = 127 time series

My aim is to cluster time series in groups which are going to show me different dynamic of sales. Before clustering I want to use Fast Fourier Transform to change time series on vectors and take into consideration amplitude etc and then use a distance algorithm and group products.

It's my first time I deal with FFT and clustering, so I would be grateful if anybody would point steps, which I have to do before/after using FFT to group dynamics of sales. I want to do all steps in R, so it would be wonderful if somebody type which procedures should I use to do all steps.

Thats how my time series looks like now:

Please note that I am relatively new to time series analysis (that's why I cannot put here my code) so any clarity you could provide in R or any package you could recommend that would accomplish this task efficiently would be appreciated.

P.S. Instead of FFT I found the code for DTW here -> link
but cannot use it on my data base and time series (suggest R to analyze new time series after 26 weeks). Can somebody explain it to me?


texts:
I have been trying to figure out how it is I have an IIR filter with large negative group delay.  It's my understanding that if a filter is causal and stable, the group delay can't be negative at all frequencies.  Is this correct?

The essence of my problem can be seen in the following MATLAB code.  This first filter gives four poles inside the unit circle and four zeros at zero.  The group delay is both positive and negative with a mean of zero.  I believe the contribution of each pole and zero adds and negates one sample of average group delay, respectively.

(1a) grpdelay(1,[1 0 0 0 .1],100000,'whole')

Now, if I create a semi-stable filter, the poles on the unit circle contribute 1/2 sample, so the average group delay is -2.  Actually it's -2 everywhere.  If the poles are outside the unit circle, there's no contribution; so, the average is -4 for the four zeros at zero.

(1b) grpdelay(1,[1 0 0 0 1],100000,'whole')
(1c) grpdelay(1,[1 0 0 0 2],100000,'whole')

OK so far.  Now if I just change the last coefficient to be complex valued, I believe this just rotates the poles, but now there is no contribution for poles inside the circle, poles on the circle subtract 1/2 sample instead of adding, and poles outside the circle subtract 1 sample.  And thus, you get negative group delay everywhere, even when the poles are inside the circle.

(2a) grpdelay(1,[1 0 0 0 .1j],100000,'whole')
(2b) grpdelay(1,[1 0 0 0 1j],100000,'whole')
(2c) grpdelay(1,[1 0 0 0 2j],100000,'whole')

I briefly glanced at the code for grpdelay but did not make a serious attempt to understand it.  Also, if you flip this around and use an FIR filter, changing the last coefficient to complex does not change the average group delay contribution by each zero.

Example (1a):

Example (2a):


texts: Estimate the Discrete Fourier Transform / Series of a Signal with Missing Samples

texts:
Assuming we have a discrete signal $ { \left\{ x \left[ n \right] \right\}}_{n = 1}^{N} $.
Which has a Discrete Fourier Transform / Series.

Now, assume I'd like to estimate its Discrete Fourier Series coefficient yet some samples of $ x \left[ n \right] $ are missing (The indices are known).

How could that be done efficiently without computing the Pseudo Inverse of the adapting Fourier Series matrix?


texts: Why is there a difference between MATLAB's FFT of a zero-padded signal compared to the non-padded case

texts:
I created a data signal and calculated the spectrum via the fft-command in MATLAB. The data-signal itself is not periodic but it is created in a way that it can be treated as one period of a periodic signal. The last sample smoothly transitions to the first one. I observe that the absolute of the spectrum strongly depends on the zero-padding before calculating the fft. If I do zero-pad the signal such that its length is a power of two, the spectrum differs severely from the one where I kept the signal unchanged. In the following pictures you can see the results.

It can be clearly seent that there is a kind of gap at the center of the channels without zero-padding. I can't explain this. I would expect the shape of both transforms to be the same. Why is that not the case here, does anyone have an idea?


texts: Fourier Transform frequency translation property in practical applications

texts:
From the material of my Signal Analysis and Processing course:

  The frequency translation property is extremely important for
  practical applications. It allows taking a signal from a portion
  of the spectrum and relocating it onto any other portion of the
  spectrum. Thanks to this property, for instance, radio
  transmission is possible, including FM radio, cell phones, TV
  transmission and so on.

I understand this property but I can't find anything on the Internet regarding its practical applications. Do you have links to any resource of this type?


texts:
I'd like to do a deconvolution of image. For example for convolution I'm using a $3\times 3$ mask with all elements $= 1$:

$$\begin{bmatrix}1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1\end{bmatrix}$$

For a pixel $X$ as follows:

$$\begin{bmatrix}X_1 & X_2 & X_3 \\ X_4 & \ \mathbf X & X_5 \\ X_6 & X_7 & X_8\end{bmatrix}$$

which is convolved by a kernel given by:

$$\begin{bmatrix}a & b & c \\ d & e & f \\ g & h & i\end{bmatrix}$$

I believe $X$ should be determined as follows

$$X=\frac{aX_1+bX_2+cX_3+dX_4+eX+fX_5+gX_6+hX_7+iX_8}{a+b+c+d+e+f+g+h+i}$$

This is a basic convolution. Now I'd like to do a deconvolution of an image which was convolved with above mask. Does anybody know how to count that?


texts:
What is meant by "spectral whitening" in DSP?

What effect does spectral whitening have when used in image processing?  (visually or otherwise...)

Where might spectral whitening be useful in audio processing or analysis?  What would a spectrally whitened audio signal sound like?


texts:
I should say which of the filters from the following image implements an ideal low-pass filter with cutoff frequency $\frac{\pi}{4}$. In the picture, Sa is a zero insertion filter, which inserts a zero sample after each input sample and Sb is a decimator, which extracts every second sample of its input. The first filter is fine and the second one does not represent the desired low-pass (that's the answer on my book).

However, I've not been able to explain why is the second filter wrong.
I know that up-sampling compresses the frequency content of a signal and that down-sampling expands it. So for the first filter, we first compress the input band by 2 and conclude that if the original limits were less than or equal to $\frac{\pi}{4}$, nothing gets filtered out (because $\frac{\pi}{4}$ becomes $\frac{\pi}{8}$ after Sa). And then, Sb recovers the original spectrum by multiplying by 2.

I analysed the second filter in a similar way, just noticing that Sb was switched with Sa thus making a multiplication by 2 first. The first Frequency Response would not filter anything from inputs with frequency lower than $\frac{\pi}{4}$ (which becomes $\frac{\pi}{2}$ after Sb). Then, after Sa, we have a division by 2 which restores the original spectrum and anything less than $\frac{\pi}{4}$ can proceed to the output.

What am I doing wrong?


texts: FIR Filter implementation without padding

texts:
Is it possible to implement FIR filtering action without padding the input and coefficients?

i.e. Let's say if the input and filter coefficients are of size 4, then the output will be 7 samples. So, while implementing, we generally add 3 more zeros to both input and filter coefficients making them equal to output size.

But, if the input and filter coefficients are of size 1024, then the output will be of 2047 samples. So, now, we need to add 1023 zeros to both input and filter coefficients. This is inefficient, right?

So, I just want to know is there any other way to implement FIR Filtering without padding?

The below code gives the idea I was talking about.

int main()
{
   int x[4],h[4],y[7]={0};

   int i,j;

   for(i=0;i<4;i++)
   {
        x[i] = i+1;
        h[i] = i+1;
   }

   for(int i=0;i<7;i++)
   {
      for(int j=0;j<4;j++)
      {
         y[i+j] = y[i+j]+x[i]*h[j]; //filtered signal of length M+N-1
      }
   }

   for(i=0;i<7;i++)
     printf("%d\n", y[i]);

}

Output: 1,4,10,20,25,24,16 - expected

But obtained results - 1,4,10,20,garbage value, garbage value, garbage value


texts: What is the meaning of half window functions?

texts: In fig. 4 on p. 4 of the paper Extracting the frequencies of the pinna spectral notches in measured head related impulse responses (J. Acoust. Soc. Am. 118, 2005 ), a half Hann window is used (Hann window for positive time, zero for negative time). What is its meaning?


texts:
I'm trying to write a simple program to perform T1 mapping of a few MRI Images.
My understanding is that all I need is to plot the pixel intensity values by inversion time then fit a curve. The output of fitting the curve is what gives me the T1 values.
This was my attempt:
# Loop through every pixel in image
for row in range(rows):
    for col in range(cols):
        x = []
        y = []
        # Add all the inversion time and pixel intesities to list
        for image in dicom_images:
            x.append(image.InversionTime)
            y.append(image.pixel_array[row][col])
        # Fit curve
        t1 = np.polyfit(x, y, 2)

Assuming what I did above is right so far I am now stuck because I don't know how to take the three coefficients produced by np.polyfit and convert them into a single pixel value. I also don't know if there are additional steps I need to take.
I found some MATLAB examples online but they were really big and hard to follow. This is my first time doing T1 mapping so I'm trying to keep it simple.
Here are some plots of the data points:


texts: direct Frequency domain FIR filtering vs Overlap-add method

texts:
I'm trying to do bandpass filtering of a EEG signal samples at 250Hz and benchmarking the following 4 methods of FIR filtering for different filter orders. The length of the signal is 15000 samples.

Time Domain Approach (convolution)
Direct frequency domain approach
Overlap-add
Overlap-save

For follow the procedure listed in the following Wiki pages for overlap-add and overlap-save
http://en.wikipedia.org/wiki/Overlap%E2%80%93add_method
http://en.wikipedia.org/wiki/Overlap%E2%80%93save_method

For frequency domain approach I do FFT for the entire length of the signal, multiply with the frequency response of filter H and then IFFT the result.

For overlap and add, I chose M(overlap) to be 1+length_of_response and L to be twice of M.

Although one would expect Overlap-add to outperform Direct Frequency domain approach, this is not what I see in my benchmark results(see attached figure).

Please help me out in understanding where I might be going wrong.


texts: When can the $\mathcal Z$-transform be inverted? When not?

texts: What are the conditions that must be satisfied to be able to invert the $\mathcal Z$-transform?


texts: Why can't Adaptive filter able to identify a system if we choose an exponential decaying input signal?

texts:
I am using the following basic block diagram for System Identification using Adaptive filter:

Here, i have used $4th$ order fir filter $h(n)$ [having weights $h=\{ 0.5 , \, 0.3 , \, 0.2 , \, 0.9 , \, 0.6 \}$ ,where $h(1)=0.5$ ] to model the unknown system for the sake of simplicity.
And i have taken the initial weights of $w(n)$ as $0$ and step size of LMS algo $\mu=0.1$.
Now, when i am using a Random Gaussian Signal as input ,then we are able to identify the system {$h(n)$};
i.e, when, $$x(n)= \frac{1}{\sqrt{2 \pi}} \exp{(\frac{-n^2}{2})}$$ with mean $\mu=0$ and $\sigma^2=1$,
we are able to identify the weights of $h(n)$ by $w(n)$ after around $50$ iterations
But when we are using a deterministic exponential signal,i.e
$$x(n)=\exp{(-n)}$$
we can't able to identify $h(n)$ by $w(n)$
When Number of iterations=$50$, $w(n)=\{ 0.0108 , \, 0.0110 , \, 0.0123 , \, 0.0196 , \, 0.0152 \}$,
When Number of iterations=$200$, $w(n)=\{ 0.0108 , \, 0.0110 , \, 0.0123 , \, 0.0196 , \, 0.0152 \}$,
and so on $\dots$
so, why can't we able to identify the system when we are using an exponential decaying signal as input?


texts: Fourier series approximation: DC component and fundamental frequency

texts:
In the linked image below, what is meant by plotting the DC component and fundamental frequency for a Fourier series approximation?

For dot point 1 does it want me to graph just the DC component basically or for the first harmonic? Basically what is meant by fundamental frequency?


texts:
I need to decode $\frac\pi4$-DQPSK packets from an SDR, preferably in C#.

I have the SDR receiver part done. Now I'm trying to learn how to take the IQ stream and get the packet info from it.

I've gone through the math, but don't understand the notation. I've tried GNU Radio, but can't get blocks added on a windows machine. I can't find a library that I can get working. I know this is rather opened ended question, but can anyone steer me in a direction to do this? I've been programming for going on 42 years now, but this is my first foray into DSP, and at my age it doesn't come easy anymore.

If it helps any, I do know a lot about the packets I'm trying to grab. They are transmitted in pi/4 DQPSK at 16kbps or 32kbps (there are two packet types), with a preamble of f337eeb637660672. The packet format is known to me, and is not encrypted.  This is for a project that will be released for free, which combined with my being retired hampers any ability to pay for professional help.


texts: My impulse response does not tend to zero

texts:
I am doing an IFFT of frequency response data achieved with Simulation tools.

When I plot my impulse response it looks wrong as the response does not tend to zero as it should. Instead there seems to be a rise in pressure towards the end of the signal.

My script looks like this:

close all;
clear all;
A = load('PvaluesAtReceiver.txt');
f = A(:,1);
Preal = A(:,2);
Pimag = A(:,3);
Prealpad = cat(1,Preal,zeros(48000-180,1));
Pimagpad = cat(1,Pimag,zeros(48000-180,1));
fpad = cat(1,f,zeros(48000-180,1));
enter preformatted text here
j = sqrt(-1);

Prealf = Prealpad(88:177);
Pimagf = Pimagpad(88:177);
ff = fpad(88:177);
Prealfpad = cat(1,Prealf,zeros(48000-90,1));
Pimagfpad = cat(1,Pimagf,zeros(48000-90,1));
ffpad = cat(1,ff,zeros(48000-90,1));
X_ff = Prealfpad + j*Pimagfpad;
X_ff = [zeros(ff(1),1); X_ff];
IRf = ifftshift(X_ff);
x_tf = ifft(IRf);
figure,plot(real(x_tf)), title('Impulse response filtered');

This produced the following plot:

Does anyone know why this doesn't tend to zero??

This is the frequency domain i want to transform.


texts: How to calculate STFT of a function for a rectangular window

texts:
How to calculate the STFT (by hand) of $$u(n)\cos(0.2\pi n)$$ for a rectangular window of a length 20, positioned at $n = 5$.

I know that to use STFT I need to divide longer signal to a shorter parts and than calculate Fourier Transform on each part. Also doing it from the definition is very long and I assume that there is a more efficient way to do it on paper.

I don't know how to start with the question so any materials on this topic are appreciated.

I came up with this idea:

First picture is just a rectangular window of length 20 and position 5.
Second picture is how I see the STFT of it. Now to calculate STFT should I provide some coefficients like mainlobe width and highest sidelobe? How I can calculate them from my data?

Are my pictures good?


texts: Manual filter or digital filter first for EEG?

texts:
I am getting conflicting advise regarding how to clean my EEG data:

1) Manually remove artefacts first and then apply digital filters

OR

2) Apply digital filters first and then manually remove artefacts

The reason given for 1) is because artefacts are more visible and avoids accidentally accepting artefacts as EEG data.

The reason given for 2) is because it avoids accidentally rejecting EEG data.

Both reasoning makes sense to me and I am a little confused as to what I should do: 1) or 2). Anyone can advise? Thank you very much.


texts:
I want to linearize a system to this form
$$\begin{bmatrix}
   \Delta\dot{x}_1\\
  \Delta\dot{x}_2\\
  \Delta\dot{x}_3
 \end{bmatrix} = A\begin{bmatrix}\Delta x_1\\ \Delta x_2\\ \Delta x_3\end{bmatrix}+B\begin{bmatrix} \Delta u\end{bmatrix}$$
    $$\begin{bmatrix}
   \Delta y
 \end{bmatrix} = C\begin{bmatrix}\Delta x_1\\ \Delta x_2 \\ \Delta x_3\end{bmatrix}+D\begin{bmatrix} \Delta u\end{bmatrix}$$

The A matrix

$$A = \begin{bmatrix} \frac{\partial \dot x_1}{\partial x_1} & \frac{\partial \dot x_1}{\partial x_2} & \frac{\partial \dot x_1}{\partial x_3}\\[6pt] \frac{\partial\dot x_2}{\partial x_1} & \frac{\partial \dot x_2}{\partial x_2} & \frac{\partial \dot x_2}{\partial x_3}\\[6pt] \frac{\partial\dot x_3}{\partial x_1} & \frac{\partial \dot x_3}{\partial x_2} & \frac{\partial \dot x_3}{\partial x_3}\end{bmatrix}_{|P}$$

$$=\begin{bmatrix} M_1L & \frac{d}{dt} \left[\frac{R_2}{M^2}\ln(x_2(t)+1)\right] & R_1\\[6pt] L_1L_2   &   \frac{d}{dt} \left[\frac{L_1}{M^2}\ln(x_2(t)+1)\right]& R_2\\[6pt] 0 &\frac1C & 0\end{bmatrix}_{|P}$$

has in its second comulmn logarithm. Pay attention only there. ($M$, $L$'s and $R$'s are some constants) If I derivate it, I get

$$=\begin{bmatrix} M_1L & \frac{R_2 \frac{dx_2(t)}{dt}}{M^2(x_2(t)+1)} & R_1\\[6pt] L_1L_2 & \frac{L_1 \frac{dx_2(t)}{dt}}{M^2(x_2(t)+1)}& R_2\\[6pt] 0 &\frac1C & 0\end{bmatrix}_{|P}$$

Am I right?

And I'm not sure what to do with this. Since my initial condition $P$ are in equilibrium and  $\frac{dx_2(t)}{dt} = 0$, should I put there the $0$? The result would look like

$$=\begin{bmatrix} M_1L & 0 & R_1\\[6pt] L_1L_2 & 0 & R_2\\[6pt] 0 &\frac1C & 0\end{bmatrix}_{|P}$$

But few weeks ago I have been solving similar situation and I didn't put $\frac{dx_2(t)}{dt}$ in the matrix at all. Lecturer told me it was right. According to that, this matrix should look like

$$=\begin{bmatrix} M_1L & \frac{R_2}{M^2(x_2(t)+1)} & R_1\\[6pt] L_1L_2 & \frac{L_1 }{M^2(x_2(t)+1)}& R_2\\[6pt] 0 &\frac1C & 0\end{bmatrix}_{|P}$$

But I'm confused. The mode $x_i(t)$ isn't consider as a function of time?

How to solve this?


texts: Help choosing coefficients for a simple FIR filter to extract DC component

texts:
I am a beginner in design of digital filters... I have just completed Coursera's Linear Circuits course!

I have a signal which is sampled at 50 Hz where the DC component has a value of 800000 and there is some noise with a standard deviation of around 200. The DC component should remain constant but it may change very slowly over time.

I want to design a simple FIR filter to obtain the DC component: would a simple weighted average algorithm be most appropriate? How do I go about choosing suitable coefficients?


texts: Time scaling and shifting of delta function

texts:
Learning signals and systems.
Solving time scaling and shifting problems.

For the  question
$$x(t) = u(2t - 1)$$
First we shift by 1 to the right side and then we do time scaling , i.e divide by 2 on the time axis.

$x(t) = δ(2t - 1) $

Can we do the same thing for the above impulse function

Is this the correct order of solving this:

Shift to right by 1
time scale by half
change the area of the delta function by multiplying it to 1/2.


texts: What is the bandwidth of a (real) sinusoidal tone, and pulse?

texts:
I would like to know how to go about calculating the bandwidth of:

A constant (real) sinusoidal tone
A (real) sinusoidal pulse.

The question is as simple as that, but I am having a hard time with the concept of what exactly the bandwidth of a constant tone should be to begin with, and from there, what the bandwidth of a pulse should be.

In the frequency domain, a constant real tone of frequency $f$ exists as two delta functions, located at $f$ and $-f$, but how does one go about calculating its bandwidth?
Furthermore, in regards to the pulse, this is rectangular function in time, and thus a sinc in frequency domain, so would not its bandwidth simply be $\frac{1}{T}$, where $T$ is the duration of the pulse?


texts: Find fundamental frequency in short time Fourier Transform

texts:
In an implementation of phase vocoder using Python, I found the author get an approximation of the fundamental frequency of STFT using the following formula:

Freq0 = B/A * abs(B/A)

where B is the fft of the first windowed frame and A is the fft of the second whose starting point is Ra(analysis hop size) behind that of B. I try my best to understand this formula but still can't figure out why it approximates the fundamental frequency. Can anyone explain this? It'd be great if detailed mathematical derivation is given. Thanks in advance.


texts: how to compute the signal passing from the low pass filter?

texts:
I am currently trying to solve this question.

Let $x[n]=\cos(\frac{\pi}{2}n)$ and $h[n]=\frac{1}{5}\text{sinc}(\frac{n}{5})$.
Compute the convolution $y[n]=x[n]∗h[n],$ and write the value of $y[5].$
Hint - Given for the question is to compute the convolution in the frequency domain first.

when I calculated the Fourier transform of $y[n]$, I got
$$Y(e^{jw}) = X(e^{jw}) \times H(e^{jw})$$
$$X(e^{jw}) = (1/2) * [\delta(\omega-\pi/2)+\delta(\omega+\pi/2)]$$
$$H(e^{jw}) =
     \begin{cases}
       \text{1,} &\quad \text{if } |\omega| <= \pi/5\\
       \text{0,} &\quad \text{otherwise} \\
     \end{cases}$$
But what I don't understand is, since the range of $X$ is beyond the range of $H$, how do we multiply the two FT's? I must be misunderstanding the concept somewhere. Could anyone help explain to me how do we find the output?


texts: What are the practically relevant differences between various image resampling methods?

texts:
Mathematica's ImageResize function supports many resampling methods.

Not being familiar with this area, beyond nearest neighbour, bilinear, biquadratic and bicubic (which are obvious from the name), I am lost.

Can you point me to some source that will explain the basic (mathematical) differences between these methods, and in particular point out the practical differences (e.g. by showing sample images where the choice of method really matters and introduces noticeable differences)?

I don't have a signal processing background, so I'd prefer a "gentle" and concise introduction :-)

I'll copy here the list of ImageResize methods for those "lazy" to click the link:

  "Nearest"  nearest neighbor resampling
  "Bilinear" bilinear interpolation
  "Biquadratic"  biquadratic spline interpolation
  "Bicubic"  bicubic spline interpolation
  "Gaussian" Gaussian resampling
  "Lanczos"  Lanczos multivariate interpolation method
  "Cosine"   cosine interpolation
  "Hamming"  raised-cosine Hamming interpolation
  "Hann" raised-cosine Hann interpolation
  "Blackman" three-term generalized raised cosine
  "Bartlett" triangular window interpolation
  "Connes"   squared Welch interpolation
  "Welch"    Welch quadratic interpolation
  "Parzen"   piecewise cubic interpolation
  "Kaiser"   zero-order modified Bessel interpolation


texts: How can I compare two pictures from the camera and tell if there are enough differences to detect movement?

texts:
I want to use my phone as a CCTV system to control a room in my house, and display an alert when something is moving.

For now, what I have succeeded in doing is to grab a picture every minute and upload them via a PHP script to my server. Now, I would like to compare the current picture and the picture 1 minute ago and detect if someone entered the room. So basically, I would need to compare differences of pixels on the picture (but also taking in account that a cloud may just say hello and changed brightness during one minute)

Does anyone have a clue on how to achieve that or some doc to read?


texts: Histogram operation for improving visual object tracking

texts:
I'm developping a people tracking with nearest neighbour data association technique and kalman filter for smoothing/predicting.

The tracker works quite good but sometimes it mismatch objects if they cross them-self.. (I think it is a common problem for a tracker).

I really believe that this tracker could be strongly improved with also histogram matching (I'm using OpenCV and I've seen an easy implementation of Bhattacharyya distance that seems quite popular in this field).

The fact is that I'd like to keep an average histogram of the tracked objects (because objects slightly changes their appearance if they rotate for example) and I'm a little bit confused on that: how to continuously adapt the average-histogram-object when new observation arrives?

observations could be different in size.. how can I handle this situation? using the same bins for each observation (even if they are 10 pixels or 100 pixels?) and sum up the histograms? using a num of bins proportional to the number of pixels?  and what about doing the average of two histograms?

I think I need all histograms with the same num of bins, and maybe use a weighted average where the weight is proportional to the num of pixels, but I'd like to ask for some more accurate advice from you!


texts: Mapping the Value of a Sample in a 2D DFT to Cycles/Pixel

texts: If I have an image and its 2-D DFT of that image, what is the mapping between the value of the DFT at (u,v), and the frequency in the spatial domain in the x and y components, in cycles/pixel? I want to be able to say "There is a peak in the DFT at (u,v), which corresponds to this frequency". I don't think the answer is u cycles/pixel, since the values shift with different image sizes, even if the image content remains the same.


texts: I tried to calculate the unilateral Z transform of x[n-2], is it right?


texts: What is the effect of noise on Shannon entropy?

texts:
Since entropy is a measure of uncertainty or randomness, intuitively we would suppose that adding noise to an image would increase its entropy since we are now more uncertain about the information of the image. I actually confirmed (partially) this assumption with the following MATLAB code where i noticed that when adding gaussian noise to a grayscale image, the entropy increases, whereas adding salt & pepper noise slightly decreases the entropy. Why is this happening and how does noise affect overall the Shannon entropy of an image?

This is my MATLAB code:

load('trees.mat');
RGB=ind2rgb(X,map);
GRAY=rgb2gray(RGB);
vecRGB=RGB(:);                                 %--vector of RGB image
vecGRAY=GRAY(:);                               %--vector of GRAYSCALE image
HVECR=entropy(vecRGB);                         %--entropy of RGB image
HVECG=entropy(vecGRAY);                        %--entropy of GRAYSCALE image
GRAY_GAUSSIAN=imnoise(GRAY,'gaussian');        %--add gaussian noise to GRAYSCALE image
GRAY_SALT=imnoise(GRAY,'salt & pepper');       %--add salt & pepper noise
vec_GRAYG=GRAY_GAUSSIAN(:);                    %--vectorize image with gaussian noise
vec_GRAYSP=GRAY_SALT(:);                       %--vectorize image with salt & pepper noise
HVECGG=entropy(vec_GRAYG);                     %--entropy with gaussian noise
HVECGSP=entropy(vec_GRAYSP);                   %--entropy with salt & pepper noise

And these are the entropies i got:

Entropy of Grayscale Image (no noise)=5.4723
Entropy of Grayscale    Image (Gaussian Noise)=7.6948
Entropy of Grayscale Image (Salt &    Pepper Noise)=5.4380


texts: Minimum sampling frequency, quantization, and bitrate calculation

texts:
An analogue sensor has a bandwidth which extends from very low frequencies up to a maximum of 14.5 kHz. Using the Sampling Theorem what is the minimum sampling rate (number of samples per second) required to convert the sensor signal into a digital representation?
If each sample is now quantized into $2048$ levels, what will be the resulting transmitted bitrate in kbps?
EDIT
The answer I get is 11 bits per sample, $2048 = 2^{11}$, is this calculation correct?


texts: What does the normalization step of the Haar wavelet transform represent?

texts:
When you perform the Haar wavelet transform, you take the sums and differences, then at each stage, you multiply the entire signal by $\small\sqrt2$.

When taking the inverse transform, you multiply the signal by $\frac{1}{\sqrt2}$ for each iteration.

What does this "normalization" really represent?


texts:
I am trying to make a panadapter for a software defined radio and I am a little bit stuck and hope someone can help me.

I have an I Q signal from my radio. (Actually I am using an IQ wav file from the internet recorded at 44800 2 channels 16bit).

If I use HDSDR the spectrum looks like this:

The display from my app looks like this:

I am doing something wrong and I do not know exactly what (I am new to this and I am not very sure I completly understand the concept behind...) I will try to describe what I am doing and hope you can help me fix it

I open the file and start reading from it in a byte buffer (1024 bytes at a time)
I convert every 2 bytes into shorts (ByteOrder.LITTLE_ENDIAN)
I apply HanningWindow to the buffer
I split the buffer in 2 buffers each for a channel (the %2==0 as
left and %2==1 right )
I create and array of Complex numbers considering the left[i] as the
real part & right[i] as imag part
I do FFT on the complex array (wavenumber table with size 512)
I do an FFT shift (from [1,2,..,n/2,n/2+1,..,n-1,n] to [n,n-1,..,n/2+1,1,2,..,n/2] )
I plot the result

Do I need to process the IQ signal before somehow (demodulate it ?)

I am reading in a wrong way the file? (I know that the file is starting with header & everything but after that the data part starts)

How can I do a pitch /amplitude correction to remove the unwanted image existing in a not ideal IQ recorded by soundcard?

Thank you,
Bogdan


texts:
I understand the reasons for zero-padding a finite signal before FFT.

However, I was reading this book (http://www.dsprelated.com/dspbooks/sasp/Zero_Phase_Zero_Padding.html) where the topic is Zero Phase Zero Padding.

Top image: windowed finite signal.
Bottom image: Zero Phase Zero Padded version of the same signal.

The merit of Zero Phase Zero Padding over "regular" zero padding (whether adding zeros before, after or at both ends of the finite signal) isn't obvious to me. What is its advantage?

EDIT:

To help anyone else who might have this question in the future, I am adding my comments here and things I tried in order to understand this more - hopefully this will save time for anyone with this question in the future.

After @matt-l's comment, I wrote a matlab script for seeing the difference in the magnitude and phase responses of a symmetric signal that was zero-padded in 3 ways (signal000000, 000signal000 and nal000000sig). The third one is Zero-Phase Zero Padding.

The imaginary part of the frequency response in case 3 is very close to 0 (on the order of magnitude of the noise floor - eps in Matlab). However, the unwrapped phase still exists - which, as Matt pointed out might be due to numerical issues.

High-res

The script I wrote is:

NSIG = 128; % Signal length
NFFT = 512; % FFT length

% Make signal
n = 0:NSIG - 1;
w1 = pi/8;
w2 = 3*pi/4; % w1 and w2 are normalized frequencies of the signal

signal = 0.5*cos(w1*n) + cos(w2*n);

% Make signal symmetric about NSIG/2
signal = [fliplr(signal), signal(2:end)]; % Make signal symmetric
NSIG = length(signal);

% Window the input
window = ones (size(signal));
% window = hanning(NSIG)';
raw_in = signal.*window;

%% A. Zero pad at the end only
NZERO = NFFT-NSIG;
in = [raw_in, zeros(1,NFFT-NSIG)];
show_spec1 = fft (in);

%% B. Zero pad on both sides
in = [zeros(1, floor(NZERO/2)), raw_in, zeros(1,ceil(NZERO/2))];
show_spec2 = fft (in);

%% C. Zero phase zero pad
in = [raw_in(floor(NSIG/2)+1:NSIG), zeros(1,NZERO), raw_in(1:floor(NSIG/2))];
show_spec3 = fft (in);

% Prepare to plot
fs = 500;
fres = fs/NFFT;     % Freq resolution
w = fres.*([0:NFFT/2, -NFFT/2+1:-1]);
figure;

%% PLOTS
titleStr = 'A. - s i g n a l 0 0 0 0 0 0 -';
h(1) = subplot (3,2,1); stem (w, abs(show_spec1),'.-'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Magnitude');
h(2) = subplot (3,2,2); plot (w, unwrap(angle(show_spec1)),'r'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Unwrapped Phase');

titleStr = 'B. - 0 0 0 s i g n a l 0 0 0 -';
h(3) = subplot (3,2,3); stem (w, abs(show_spec2),'.-'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Magnitude');
h(4) = subplot (3,2,4); plot (w, unwrap(angle(show_spec2)),'r'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Unwrapped Phase');

titleStr = 'C. - n a l 0 0 0 0 0 0 s i g -';
h(5) = subplot (3,2,5); stem (w, abs(show_spec3),'.-'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Magnitude');
h(6) = subplot (3,2,6); plot (w, unwrap(angle(show_spec3)),'r'); axis tight; title (titleStr); grid on;
xlabel('Frequency (Hz)'); ylabel('Unwrapped Phase');


texts: relationship of frequency with range in Underwater acoustic channel

texts: I am interested to know what range of frequencies can be used/ or are effective for wireless communication, within the underwater acoustic channel. I do know a table exists from \cite{akyildiz2004challenges, akyildiz2005underwater, kheirabadi2013greedy}, which is reproduced below. I want to add another column in this Table to add the range of frequencies that would be effective for the given range of propagation distances. But cannot get an idea of how to get it done.


texts: I am trying to find a specific block of data that should hold valid audio data that I want to process. For instance, I record at 44.1kHz, but I only transmit a .5s of content. Since the buffer I store the recordings in contains both the recording as well as background noise, is there a way to isolate the sample I desire? Is this a simple filtering problem? I know the range of frequencies I need in the sample, so would a bandpass filter solve this issue? My ultimate goal is to run the time series data through an FFT to see which of the desired frequencies were transmitted. Note that this is being done in software. Thanks for any help!


texts:
i have database for 50 person ,it is 250 images (5 images/person) , i need to perform one to many matching between images for recognition persons.

so , i divided the image as 2 parts (one for training and one for testing)

for example if we say that the first folder contains 100 images(2 images /person),and the second contains 150 images (3 images /person)

is this division correct?
 if it is correct , how can i matching each image in folder1 with all image in folder 2(which will contain 3 similar image) and obtain the better result after (that give max matching percentage)from folder 2 and finally returning the name of this image ?

i need to perform this in matlab

please , i need many advices and consultations with any one has good knowledge in this field.

this is an example of my images:

regards


texts: Why does this constant-Q transform "lose octaves"?

texts:
In this paper (Schörkhuber, Klapuri, Sontacchi), it's shown that the constant Q-transform evaluated "octave-wise" produces the following looking graph:

I don't understand why the octaves are "disappearing" in some time intervals. Isn't this then a lossy" process somehow?

In the paper it says that this is achieved because

  the atom hop sizes $H_k$ are set according to the highest frequency bin within each octave, that is, for each octave down $H_k$ is multiplied by 2


texts: Calculate bandwidth of a signal for Nyquist–Shannon sampling theorem

texts:
I have to calculate the minimum sampling frequency for the Nyquist–Shannon sampling theorem which is Fc > 2*B, where B is the signal bandwidth.
I have this signals: $\text{sinc}^5(t/2 - 4)$ and $\text{sinc}^3(3 - 2t)$

How can I calculate the signal bandwidth to obtain the minimum sampling frequency for the Nyquist–Shannon sampling theorem?


texts: Is power density invariant to Fourier Transform? Does it hold through derivation?

texts:
Assumptions

We have a finite discrete measurement.

Introduction

I had some troubles determining a power spectral densities needed for Wiener deconvolution. Looking through the continuous equations of FT applied to power and spectral power it seemed like they should be basically equal. Then I got distributions involved, resulting in following proof.

Let $<P|\varphi>$ be a distribution (linear continuous functional on set of testing functions $\varphi$).  Define $P$ as $$P(\varphi)=\int_{-\infty}^{\infty}{\varphi^{2}(x)dx}$$
Let us have a Fourier transform of a distribution (can prove this one), $$< F(P) |\varphi>=<P|F(\varphi)>$$
Then let us demonstrate the equality:$$<F(P)|F^{-1}(\varphi)>=<P|F(F^{-1}(\varphi))>=<P|\varphi>$$
Plugging in the definition of the distribution, we get:$$P(\varphi)=\int_{-\infty}^{\infty}{\varphi^{2}(x)dx}=\int_{-\infty}^{\infty}{F(\varphi)F^{-1}(\varphi)dx}=\int_{-\infty}^{\infty}{F(\varphi)F^{*}(\varphi)dx}=\int_{-\infty}^{\infty}{{|F(\varphi)|}^2dx}$$ Thus proving the invariance of total power included in the signal.
This holds on through simulations (ratio between total power and total spectral power converges to 1).

SNR ratio

Wiener deconvolution uses a ratio of power densities of a signal and some noise. For a fixed range over which we assume the density, explicitly the range on which the signal is non-zero, we should be able to determine explicit ratio needed. Using all the assumptions, I was able to contradict this approach via derivatives.

Let us have a derivation of a distribution (also provable):$$<P'|\varphi>=-<P|\varphi>$$ Therefore the negation should cancel itself in the desired ratio, however it does not (using the same simulation, it simply varies).

Question

Somewhere I had to make a mistake. I suppose it happens in continuous theory when I subconsciously multiply by Dirac's comb, or when I assume the same length in the density formula $$\int_{-a}^{a}{\frac{f(x)dx}{2a}}$$

Plea

I hope I just need another pair of eyes, because if a theorem holds, a result is derived and observed, a theorem is used again and it contradicts itself, there must be something wrong...


texts: What is the relationship of S-parameters to frequency response?

texts:
I have some S-parameters. I'm reading about frequency response and the impulse response and step response.
If I have a set of S-parameters, then is the frequency response for a port (say output port 2, input port 1) simply $S_{21}$ for all the frequencies in my S-parameter file?


texts: Why are Gold codes and Kasami codes used instead of pure m-sequences?

texts:
I'm trying to understand why Gold codes and Kasami codes are used instead of pure m-sequences, in direct-sequence spread-spectrum (DSSS) communication systems, to prevent interference between multiple transmitters.

If I understand correctly, Gold codes are defined as the XOR between two m-sequences with different polynomials of the same degree (e.g. one LFSRs based on $p_1(x)$, and another LFSR based on $p_2(x)$ both of degree 20; the LFSR outputs are then XOR'd together), and within one system, multiple transmitters use the same pairs of polynomials but with a different time shift. (Kasami is also an XOR so I guess it's similar.)

How is this better than, say, having a system where all transmitters use m-sequences but each of them uses a different polynomial? (In both cases, the autocorrelation values are higher than the cross-correlation values)


texts: is it possible to introduce fractional delay in waveform in GNURadio?

texts: I know this is a strange question as its hard to do in digital domain. I know it is a straight forward task if the signal is defined within N number of samples where I can go to Fourier transform multiplied by complex sinusoidal with inverse period of that fractional delay and then come back to time domain through IFFT. Is it possible to do in GNURadio ?


texts: Filtering out unwanted square wave (Radio: DCS/DPL signal)

texts:
I've need to remove a DPL/DCS signal from a recording. Roughly speaking, a DCS signal is a low frequency (67.15Hz) square wave sent at the same time, over the same carrier, as speech. Because it is a square wave, it has many strong harmonics that overlap with speech. Obviously, the speech must be preserved.

On the surface, this seems like a problem that might be solved the same way as removing 60Hz power-line noise: lots of notch filters. However, power-line noise tends to be weaker and comes from a source that is much closer to a sine wave.

So, my question is: is there a better way to do this?

This link contains more info about DCS: http://onfreq.com/syntorx/dcs.html
It mentions "Since DCS creates audio harmonics well above 300 Hz (i.e. into the audible portion of the band from 300 to 3000 Hz), radios must have good filters to remove the unwanted DCS noise." Ha!

EDIT: DCS is not a "square wave". I think the correct way to explain it is to say it's a NRZ binary code clocked at 67.15 Hz.


texts: I'm trying to turn audio into haptic feedback, how do I get numerical data from audio?

texts:
How would I sample a song to get numerical data?
I was thinking there might be some kind of program that could do this or maybe I could interpret an audio file itself without processing it in any way.
I'm going to start with 3 haptic feedback vibration motors, one on each arm and one on the chest.
I would probably need to break that data up into different parts such as stereo left/right and base.
If I can't get multiple channels perhaps I could play different frequencies on different motors.
Once I have the data, I would have to figure out how to interpret it. I think I can figure this part out. My plan is to then take that data and copy it to a Raspberry Pi and "play" it on the vibration motors.
I don't know how fast I could get the motors to accurately go. I've played around with them a little. Maybe 1-second intervals so I may have to skip some beats in faster songs.


texts: What is importance of the conjugate reciprocal and none reciprocal zeros in FIR Type 2 Filters?

texts:
In DSP class we have this slide:

So i like to know why this classification is so imoortant?

I have seen this post :

  Conjugate reciprocal pairs of zeros and poles in ...

  Assuming the impulse response h[n]h[n] of an FIR filter is real for
  all nn,

  Why are zeros and poles in FIR design found in reciprocal and
  conjugate pairs?Is the assumption necessary for this phenomenon to
  take place?

But i dont get my answer?( mybe because of being newbi in this field)


texts: Python vs Matlab? Which one better for image processing?

texts: I want to learn practical image processing from scratch/basics in easy way. Which software will be helpful in this regard? Matlab or python?


texts: Is sparsity induced penalty in source separation "Entrywise matrix norms"?

texts:
I am reading this paper where they introduce norm penalties for source separation. In table 1, the $\log/ l_1$ type is $\sum_{g} log(\epsilon +  \lVert H_{g} \rVert_1)$. I wonder this $\lVert H_{g} \rVert_1$ refer to which definition of $l1$ norm? In wikipedia, it has 2 definitions for this symbol.

"Entrywise" matrix norms: $\Vert A \Vert_p = \Vert \mathrm{vec}(A) \Vert_p = \left( \sum_{i=1}^m \sum_{j=1}^n |a_{ij}|^p \right)^{1/p}$ (wiki) or
Matrix norms induced by vector norms: $ \|A\|_1 = \max_{1 \leq j \leq n} \sum_{i=1}^m | a_{ij} |$ (wiki)

which is simply the maximum absolute column sum of the matrix.

One thing is below the table, the description said ""all norms are elementwise over matrix entries". So I doubt it is the first definition but not fully understand. Very appreciate if you may add explanation why it is that type of norm.


texts: Is there a function in MATLAB similar to $\tt freqz$ but for $H(\omega)$?

texts: I wanted to know if there's a MATLAB function that equates the frequency response, $H(\omega)$ of a system. Something similar to freqz like a freqw.


texts: Amplification of a speech signal in the recorded audio

texts:
I have recorded an audio.
I dont know how it happened that only one sided speech is recorded and the other speech is recorded with a very low sound. Is there any solution to amplify the other side signal.  I have tried to amplify with AVS software but its also ampliying the noise in the audio.
any help would be much appreciated.


texts:
I am trying to implement a binary offset carrier (BOC) direct sequence spread spectrum scheme in MATLAB. In that scheme, a code $c(t)$ (a binary vector which repeats after a certain number of zeros and ones) should be calculated by an other square (sub) carrier of the following form sign(sin(1*pi*fs*t) where fs is the subcarrier frequency and finally should be calculated by a carrier frequency having the usual form of a cosine.

I have generated the code needed (a vector of 1023 elements). These 1023 chips (vector’s elements) should be repeated every 1 ms in order for the chipping frequency to be achieved (10Mchips/s). After, that vector (code) should be multiplied by a square subcarrier of 1023 Mhz and finally by the ‘actual’ carrier (a co sinusoidal signal) of 4092 Mhz. All these should be demonstrated via plots.

My first problem is that in order to demonstrate the spectrum of the codes, spectral nulls should be represented at frequencies multiples of 1.023 Mhz. The only way that I could somehow think how to do this is to use massive vectors where are not really supported by my laptop's memory (RAM).

I say "somehow", because the sampling rate should change (for example wherever I have one ‘1’ I should add some more and the same applies to zeros). This procedure ruins the whole purpose of the code because the code has some particular autocorrelation properties which I don’t know if will continue to exist after.

Moreover, if that should reach the 10Mchips how will this could be done by using the particular vector. While I understand that it is easily done by using a random vector of 10*10^6 elements. Should I repeat the 1023 chips 1000 times?


texts:
Hi I am trying to derive the autocorrelation matrix and I am unsure about how exactly it works. I can see that the $4\times 1$ matrices result in the Hermitiain and Toeplitz matrix? Surely the only non-zero values in each of these is at $n = 0$ which is the variance?

Any help would be greatly appreciated!


texts: Applying a filter on an input sequence in matlab

texts:
Given a input sequence $x = {x_1, ... x_n}$ and a filter with l elements $h = [h_1, ..., h_l]$ with l < n. We want to filter the input sequence with the specified filter.

My first thought was to compute the convolution between $x$ and $h$:

$y = \sum_{k=0}^{l-1} x(k) h(n - k)$

For this to work I'll have to pad x with zeros.

x = [ zeros(1, l) x zeros(1, l - 1)];

And then

for i = 1:length(x) - 2
     y(k) = h * x(k:k + l)';
 end

Is this the right approach for this problem ?


texts: Octave: Finding the mean fundamental frequency given an audio signal

texts:
I have an audio .wav file of someone's speech of which I want to find the mean fundamental frequency. Till now, all I have managed to do is load the data into a vector and (roughly) filter the silence out. I referred to a MATLAB example of finding the fundamental frequency, but it lies way out of the range of the mean fundamental frequencies which I have in a dataset I am referring to. The code on MATLAB is so-

function F0 = fundafreq(data)
  % Copied from matlab site
  F0 = 0;
  Fs = 48000;  % My sampling frequency
  [b0,a0] = butter(2,325/(Fs/2));
  xin = abs(data);
  xin = filter(b0,a0,xin);
  xin = xin-mean(xin);
  x2 = zeros(length(xin),1);
  x2(1:length(data)-1) = xin(2:length(data));
  zc = length(find((xin>0 & x2<0) | (xin<0 & x2>0)));
  F0 = 0.5*Fs*zc/length(data);
end

This gives me the fundamental frequency of around 40 Hz whereas my data set has mean fundamental frequencies not more than 1. I am only a beginner when it comes to digital signal processing, so I am not completely aware whether this code is giving me the answer or whether I need further calculations.

Also, it would be preferable if I could calculate this without needing the sampling frequency.


texts:
I am studying proakis book digital signal processing using matlab 3rd ed

but i am bit confused about calculation of value complex multiplier W

fig is attached. i am not able to understand how the values of W are found/calculated in this fig in red enclosure

how/why we know that W base 4 and raised to the power 1 is equal to -j


texts: Improving an auto-correlation based guitar pitch detector

texts:
I've seen many questions on this forum regarding pitch detection for musical instruments (commonly guitar), and spent a while reading through the answers to create a basic implementation of auto-correlation to make an Android guitar tuner.

This is the algorithm I'm using (implemented in Java on an Android phone):
1) Record a short[4096] array of raw audio data from an Android phone's microphone at 48kHz sample rate
2) Apply a Hanning window to the raw audio data
3) Zero-pad the result to double the length (8192)
4) Apply auto-correlation with FFTs
5) Normalize the auto-correlation result
6) Get the periodicity from the peak bin indexes

My problem is that with an actual guitar it is not robust (around 50% accurate at best), and I don't know how to filter noise either (without any loud noise, just ambient white noise, it outputs garbage frequencies).

It fares much better when I whistle at it, or play a generated sine tone from a computer, but that's expected.

In my search for ways to improve my implementation, many answers (on this forum and others) usually point to using a better algorithm like YIN, but I think its more likely that I've made a mistake in my implementation.

Does anything seem obviously incorrect from the algorithm I posted? Are there tweaks I can apply?

1) Increase the Android microphone sensitivty
2) Band-pass the raw data before the algorithm
3) Record more samples
4) ??? And does anybody have any idea how I could filter out noise by using the auto-correlation result?

Thanks in advance.


texts:
Can you suggest any good introductory resources on ridgelets and/or curvelets for image processing purposes?

I'm considering buying Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity (more). Is there anyone who can recommend for or against it?


texts:
I have searched the web alot but I couldn't find a page that has taught directional filtering comprehensively. Could you please learn me the subject and list the kernels for all directional filters?
Or introduce me a page that I can learn the subject from?


texts: On the signal recovery in OFDM after estimation

texts:
Given a modulated signal X which is transmitted based on OFDM with N = 256 sub-carriers through channel h = [0.7 0.3 0 0 0 0.4] of length L = 6; following the steps : [X --> ifft(X,N) --> adding CP --> ... ext] and then we have Y representing the received signal representing the output of FFT at receiver as below :

X = The modulated signal of length 256;

h = The channel of length 6;

Y = the received signal after taking the FFT of length 256 too.

Assuming the channel is estimated at receiver perfectly and h_est' = [0.7 0.3 0 0 0 0.4];

I's asking how can we get the sign X_est back ? Is it by X_est = Y. / fft(h_est, 256);    however when I try that way in Matlab, it doesn't give the estimated transmitted signal !

her is the code as reference:

clear; close all; clc

n_subc = 2^8;                          % # of subcarriers
n_ofdm_sym = 1;                        % # of OFDM symbols
n_data_sym = n_subc*n_ofdm_sym;         % # of data symbols to transmit
CP = n_subc/4;                          % CP length

M = 4;                                  %Modulation order setting

%Set channel
h = [0.7 0.5 0 0 0 0.2];

iter = 10^2;    %% Iteration
FFT = dftmtx(n_subc);           %% FFT matrix

SNR = [0:2:40];

for i = 1:length(SNR)
    i
    for j = 1:iter                                % Monte Carlo iterations

    d_symb = randi([0 3],1,n_data_sym);            %generate the data to transmit
    x_mod=qammod(d_symb,M) / sqrt(2);           %modulate

    x_mod_s2p = reshape(x_mod,n_subc,[]);        % serial to parallel

    x_ifft = FFT'*x_mod_s2p;             %Ifft conversion

     x_cp = [x_ifft(end-CP+1:end,:); x_ifft];     % add CP

    x_ifft_p2s = reshape(x_cp,1,[]);             % parallel to serial

    % --------------------> Channel mixing and adding noise
    y = conv(x_ifft_p2s,h,'same');                          % Linear convolution
   % y = awgn(y,SNR(i),'measured');                          % add noise
    % --------------------> Receiver

    y_s2p = reshape(y,n_subc+CP,[]);             % serial to parallel conversion

    y_cp = y_s2p(CP+1:end,:);                    % Removing CP

    y_fft = FFT*y_cp;                            % FFT

    y_p2s = reshape(y_fft,1,[]);                 % Parallel to serial

    X_est = y_p2s ./ fft(h,n_subc);
    y_demod= qamdemod(X_est,M);

    [~,ber(j)] = biterr(d_symb,y_demod);
    end
    BER(i)=mean(ber);
end

figure
semilogy(SNR,BER)


texts: I am working on a color transformation project using video.Can any one please tell me how to read video in octave.I have see aviread can be used to video frame by frame but i want to read the entire video file.


texts:
How many stable and causal systems with the same magnitude response are there?
I know this relates to an all pass system for two rational transfer functions but am not sure about the specifics of this


texts: Basic question regarding third octave filters

texts:
I would like to use third octave filters on a time signal.
Do I have to apply all third octave filters to a single frequency or do I only need to apply the one that has its center frequency the closest to the frequency that I am looking at?

I would like to use a simple example to make my quesion clear:
Let's say I have a basic sin signal: signal = sin(2*pi*f*t) with f = 57Hz.
57Hz is almost right between the center frequencies of 50Hz and 63Hz, but the closest would be 63Hz. Do I apply all filters or only the one with a center frequency of 63Hz?

I am pretty sure that the answer is that I will always have to apply all filters in theory. Of course filters with much higher center frequencies will result in insignificant amplitudes...


texts: Degraded image: Wiener Filter vs Wiener deconvolution

texts:
I am studying methods for recovering images (that has been blurred and noised).

Can someone explain the difference between:

http://en.wikipedia.org/wiki/Wiener_deconvolution
http://en.wikipedia.org/wiki/Wiener_filter

I can understand that wiener deconvolution is used for Deblurring while wiener filter is used for Denoising.

Is that right? Are there any difference between the two?


texts:
I have certain doubts in the fundamentals of AC and DC signals.

Most of the time, we deal with DC signals in the form of batteries, whose waveform is a constant signal with peak as 1.5v

However, i came across some resources, that told certain other signals are DC as well. For your reference, i have enclosed the resources and the respective waveforms here.

Figure 1

Figure 2, Please ignore the "bidirectional" word on the image.

So, are these signals AC, or DC? Both the sources have mentioned that they are DC signals (well, the first one is apparently "DC Pulsating" and the second one is a "DC time-varying". This confused me, and i wanted some better clarification and knowledge on these. From what i observe,

1) Any signal that is above or below the time axis (i.e the same polarity all the time) can be called DC?

2) If i have an AC signal, say a sinusoid, and by appropriately DC shifting it above / below the time axis, can i make it a DC signal, if the above point is true?

I would seek your expertise and your insights on these images as well, and tell me why these signals are DC.

Thank you.

Reference 1: https://www.elprocus.com/variable-power-supply-digital-control/

Reference 2: https://www.youtube.com/watch?v=fXBtTI107fQ, time: 2:31


texts: What is meant by "Cubic Difference Frequency"?

texts:
The cubic difference frequency (or 3rd degree intermodulation product) is given by:

$$D_3 = - 10\log\left(\frac{P_2\left(2 f_{1,2} - f_{2,1}\right)}{P_2\left(f_1,f_2\right)}\right)$$

where $P_2$ is the output power. As far as I understand, the 3rd degree IM products are at frequencies $2f_1 - f_2$ and $2f_2-f_1$. So the formula for the IM product of $2f_1 - f_2$ should be:

$$D_3 = - 10\log\left(\frac{P_2\left(2 f_{1} - f_{2}\right)}{P_2\left(f_1\right)}\right)$$

When the power is given in an logarithmic scale, this would simply be the difference between the two power values.

But what does this actually tell me?
Is this the power of my 3rd degree intermodulation product?
(if it is, $D_3$ should be equal for both frequencies)
And is this formula only valid for an input signal of two different frequencies?

According to a paper D3 can also be expressed in terms of gain compression c and change in AM-PM characteristic kp due to only one carrier:
$$D_3 = -10\log(\frac{c_1}{2}^2 + k_{p1}^2)$$

The gain compression is given by $c = 1 - \frac{\Delta P_2 / P_2}{\Delta P_1 / P1}$
, where $\Delta P_1$ is a constant change of input power and $\Delta P_2$ is a not constant change of output power. I am not sure however, if $P_1$ or $P_2$ refers to the input/output power before the change or after.


texts: Autocorrelation sequence in terms of Fourier transform of the underlying signal

texts:
Let $x(n)$ be a sequence of length $N$, which is zero outside the interval $(0,N-1)$. Let $X(k), k=0,1,\cdots,N-1$ be the FFT coefficients of $x(n)$, that is, $X(k)=\sum_{n=0}^{N-1}x(n) \exp\left( -\frac{j2\pi k n}{N}\right)$. How to relate the autocorrelation sequence  $y(l)=\sum_{n=0}^{N-1-l}x(n)x^{*}(n+l)$ in terms of $X(k)$?

My try:

\begin{equation}
\begin{split}
y(l)& =\sum_{n=0}^{N-1-l}x(n)x^{*}(n+l)\\
& =\sum_{n=0}^{N-1-l}x(n) \left[  \sum_{k=0}^{N-1}X^*(k) \exp\left( -\frac{j2\pi (k+l) n}{N}\right) \right]\\
& = \sum_{k=0}^{N-1}X^*(k) \left[\sum_{n=0}^{N-1-l}x(n)    \exp\left( -\frac{j2\pi (k+l) n}{N}\right) \right]\\
\end{split}
\end{equation}

How to proceed from here?


texts: Linear response function for a system with derivative: $U=L \frac{d I}{dt}$, expressing $U=f(I)$

texts:
I have a super basic questions. I am a not really into signal processing (more about physics), but I would like to understand an aspect of linear response function (I think the question fits for this forum).

From my basic understanding, any linear, time invariant, causal system can relate the output $S(t)$ to the input $E(t)$ under the following relations:

$$\begin{equation}S(t)=S(-\infty) + \int_{-\infty}^{+\infty} \chi(t-t') E(t') dt'\label{eq_1} \end{equation}$$

The causality holds when imposing $\chi(u<0)=0$.

Let's consider a very simple case: voltage around an inductance. I have the law: $U=L \frac{d I}{d t}$

If I express $I$ as a function of $U$, I can write down:

$$I(t)=I(-\infty) + \int_{-\infty}^{+\infty} \frac{U(t')}{L} dt'$$

My response function will simply be $\chi(u)=\frac{\Theta(u)}{L}$ (where the heaviside is for causality).

But the system in which $I$ is the input and $U$ the output is also linear. Thus I would expect to be able to express:

$$U(t)=U(-\infty) + \int_{-\infty}^{+\infty} f(t-t') I(t')$$

However because of the derivative on $I$ in the law, I don't find how it is possible.

Am I wrong in my initial statement and \ref{eq_1} ?


texts:
The complex exponential form of cosine

$$\cos(k \omega t)  = \tfrac{1}{2} e^{i k \omega t} + \tfrac{1}{2} e^{-i k \omega t}$$

The trigonometric spectrum of $\cos(k \omega t)$ is single amplitude of the cosine function at a single frequency of $k$ on the real axis which is using the basis function of cosine?

The complex exponential spectrum of $\cos(k \omega t)$ has two amplitudes at 1/2, one at $k$ and $-k$.

I am confused what this x-axis is representing - amplitudes of what? Whats the basis function? I'm guessing amplitudes of exponential but these complex exponential have a cosine and sine how can they be on one axis if they are made of cosines and sines and cosine and sine are orthogonal and orthogonal are on different axis.

Appreciate your advice, thanks.


texts: How to generate the filter bank for a wavelet adapted to a pattern?

texts:
Some textbooks like Daubechies (1991) and Strang & Nguyen (1996) give algorithms to construct wavelets with specific properties. But, I am trying to construct a filter bank for a wavelet that is as close as possible to a given pattern (vector). Seems like this approach is called matched wavelets.

Do you know if there is some code available for that?


texts: What sonic artifacts are caused by phase vocoders?

texts:
When using a phase vocoder to alter the spectral envelope of a signal, the input signal is split up into many (often hundreds of) bands. With many phase vocoders, the output levels these individual bands are often controllable by the user.

What artifacts are caused by filtering these individual bands in a phase vocoder?
If and why, phase vocoders are better suited (and cause less artifacts) than FIR filters (or other methods) for surgical changes to such narrow frequency bands.  For example, are phase vocoders going to cause the same, or less, pass band ´ringing´ than using an FIR filter for the same task?

I understand phase vocoders only up to the point that they use an inverse DFT to resynthesise the sound, though Im not sure exactly how the individual bands are resynthesized, so I dont know what the problems/artifacts will be.


texts: Sum of Squared Intensity Difference for block matching with ORB features, c++

texts:
I am trying to implement feature matching algorithm between image reference and image current.

The pipeline is as follows:

Gray scale Images (intensities = [0-255]) are imported.
Extract ORB keypoints in both reference and current images:

Ptr<FeatureDetector> detector = ORB::create();

Obtain blocks surrounding all keyPoints within specified window_size in both reference and current images. NOTE: If the blocks exceed the image upper or lower bounds, the block intensities would obviously be zero (Example: KP_A in image reference).

4.Instead of using existing feature matching algorithms in opencv, I am trying to utilize sum of squared of intensity differences (SSD) in the blocks acquired across ORB keypoints in reference and current images.

Extracting correct features demands implementing crossCheckedMatching() to ensure features are chosen correctly.
Visualize matched features.

Complete code: here.

Problem 1: My final result looks something like this:

Whereas, I have read in few lecture notes and articles, SSD and block matching incorporated with cross checked matching method should provide decent matches. Does these matched features look any good in someone's eyes or awful after cross check matching procedure?

Problem 2: I do not know how to find the correct threshold for evaluating SSD values between two images (Empirical tuning vs. Existing logical approaches?).

Final goal, adapted from the textbook:

use these matched features in RANSAC to filter outliers.
Find Homography (Theory says, 4 pairs of matched points needed!)
Camera Pose Estimation
3D reconstruction of environment

Apologies for such a long question, but I tried to squeeze it as much as I could.


texts: Can the FFT reduce its computation load for specific sparse data patterns?

texts:
I am working on a signal processing application, which involves an IFT.
My background is CS and not signal processing maths.
Consider an N-FFT configured with 256 "bins" (N = 256) for processing a
100Hz bandwidth. The application is in the a-priori position of knowing at
a given time where the input data can be sparse, in terms of only a subset
of the total bandwidth being active.
Ex: at time T, there is a signal in a band B of [26,40]Hz, while at [0,25)
and (40,100] Hz the signal is zero.
My questions are whether :

in such cases the 256-FFT could be configured to operate only
as a 16-FFT (as B = 40 - 26 = 15, 2^2 < 15 <= 2^4) , and potentially avoid
doing reams of computation as a 256-FFT where the add/multiply would be on
zero values (and thus redundant) .

if so, professional impls of the FFT already do this (a quick scan of
the input vector etc to determine whether such sparsity exists etc before
beginning the computation proper) .


texts: Get the lag of a signal with respect to the other with cross correlation

texts:
I have 2 vector of length 20 each. The first 15 entries of $a$ are exactly the same as the last 15 entries of $b$. There is random noise in the 5 last entries of $a$ as well as the 5 first entries of $b$. I would like to use cross correlation to detect the lag of 5 between those 2 signals and wrote the following code:

%Matlab Code
clear all; close all;
iteration = 128;
suma = 0;
for it = 1:iteration
  offset = 5;
  range = 100;
  vecLength = 20;
  a = randi(range,1,vecLength);
  b = zeros(1,vecLength);
  b(offset:size(b,2)) = a(1:size(b,2)-offset+1);
  noise = randi(range,1,offset-1);
  b(1:offset-1) = noise;
  c = xcorr(a,b);
  %[c,lags] = xcorr(a,b);
  [argvalue, argmax] = max(c);
  argmax = argmax - length(a);
  suma = suma + argmax;
end
lag = suma/iteration

However, I get an average of around -2.0 for the value of lag. Here are the result of 7 runs:

-1.7344
-2.1875
-2.1719
-2.1641
-1.9531
-2.3750
-2.1172

How can I get the back the value of offset = 5 instead?


texts:
I have tried the following methods to test the orthogonality on Fixed point FFT implementation.

FFT:

Y = A*x, where A = FFT basis and x = input, Y= FFT output

If the input (x) is set to A-1(A inverse), then the o/p will be an identity matrix.  To get A-1(A inverse), first set the input (x) an identity matrix. The output is nothing but the basis vector. Take the transpose and set this as next input, which should give an identity matrix as output.

16bit fixed point implementation, an input value of 1 is set as 0x7FFF, but I was not able to get the orthogonality test by setting this value. Please help me in doing orthogonality test on fixed point FFT implementations.

-ben


texts: how to detect corner in a digital image using MATLAB?

texts: I am designing an algorithm for image processing. I am using MATLAB for the processing of the images(video). How I can detect the position of the corners in the image?


texts: Does waveform-based ASR with lowpass filtering preprocessing still require the assumption of FFT?

texts:
I read some recent papers about time-domain speech recognition, and wondering if it steers away from the stationary assumption of FFT. However, some preprocessing operations like pre-emphasis and low-pass filtering still employ FFT. Does it mean waveform-based ASR with these preprocessing still has the FFT assumption of stationarity?

Many thanks!


texts: Choice of Gaussian kernel parameters when lowpass filtering before image resampling?

texts:
I need to decimate a signal by a factor of q.

More specifically my signal is a 3D "image": $\ I(x_i,y_j,z_k)$,
which I need to downsample by a factor of two in the z direction.

I want to do lowpass filtering before decimation by convolving with a Gaussian kernel of size n.

I create my Gaussian kernels 2 standard deviations below and above 0 since this accounts for 95 % of the distribution.

I am looking for a rule of thumb that tells me how large n should be.

Am I right in thinking that the Gaussian should filter out all frequencies above fN/q,
where fN: Nyquist frequency of original signal?

I know that the Fourier of a Gaussian with standard deviation $\ \sigma $ is another Gaussian with standard deviation  $\ \sigma^* =1/\sigma $. I am not sure how strict my lowpass filter should be. Should the cut frequency, fc, be at 2 or 3 standard deviations?

What is then the equation for the cut frequency, fc, of a Gaussian kernel with size n: fc(n)=?

Below is the frequency response of some Gaussian kernels calculated in Matlab:

My actual problem involves q = 2, and from this figure I see that n=5 should work nicely.
It would have been good to have a rule of thumb do, so I did not have to do this for each q I encounter.


texts:
I am trying to code audio signal processing software in java.

Purpose: It will use audio signal processing library to filter out noise from the raw PCM based input (.wav) of heartbeat recording.

Steps involved

Select input .wav file.
Apply Low Pass filter at 33 Hz.
Apply High Pass Filter at 2 Hz.
Apply Amplification with Max Gain.
Store back it as .wav file

There are many C++ library which will cover most of the requirement. But looking for Java for this particular project. Need input from Audio Signal Processing engineers for their input if these steps will be sufficient to able to hear the heartbeats. Also please let me know which order of filtering and which window will be most effective in this scenario.

Input is taken from

http://www.lehman.edu/faculty/hoffmann/itc/techteach/audio/sound2.html


texts: How to sync audio and video signals in MATLAB?

texts:
I have a video clip that has a frame rate of 25 fps (MJPEG compression) and an audio clip with a sample rate of 48kHz. Both of them play for the same duration - 27 seconds.

Is there any standard way to sync the two signals of different sampling rates if I wish to write them to file?

If this is not doable cleanly in MATLAB, I am happy with a some software tool on Linux or MacOS that can do this syncing.

E.g. the VideoFileWriter object can take both an audio sample and a video frame at a single step, however I have about 1917 audio samples corresponding to a single video frame.

Some details:

I have read in the audio signal using [signal, fs] = audioread(filename).
signal is a 1336321 x 1 array.

All the 697 frames in the video are saved as a cell array of RGB images.


texts: Cross-correlation of two signals using matlab

texts:
I am trying to measure the similarity between two signals and I am using cross-correlation to achieve this.

Following the matlab example, I seem to have accomplished that, however, I do not understand the plot that I am getting. In the matlab example the explanation that is given is: The first subplot indicates that the signal and template 1 are less correlated while the high peak in the second subplot indicates that signal is present in the second template.

My question is how does the peak indicate that the signals have more in common? Shouldn't the two signals be more correlated when the distance between them is lower?

Moreover, when I tried it using my two signals I received the following image:  which looks more than strange compared to the image from the example, as I have multiple peaks and a relatively long straight line.

EDIT

I added two example signals that yield a relatively same result. As suggested from the comments, I tried without upsampling the signal and making them equal in length, however, the results weren't appealing as well. Furthermore, if the peak is, for example, at the 10s, does this indicate that the signals are simillar at the tenth second? What about if the peaks are in the negative time zone (t<0s)?


texts:
I am trying to measure phase difference between two received waves which are recorded by two separate PC's (microphones). The two PC's starts recording at the same time, then I make a sine wave sound with frequency of $1000\textrm{ Hz}$ using a sound speaker which is recorded by the two microphones.

For the purpose of measuring the phase difference between the two recordings, I do a cross-correlation between the two waves. Theoretically ... the distance between zero and the index of the peak value of the output graph (after cross-correlation) is the phase difference between the two waves.

For some reason I am getting a graph with more than one peak. I think this is because of the noise.

Does someone know if I can make a normalized cross-correlation to sound waves, does this solve my problem  ?
If yes, how can I do this ?
If no, any suggestions to solve this problem ?

The waves and the correlation after plotting in MATLAB:

I know the two waves are not perfectly the same, they represent a sine wave with $1000\textrm{ Hz}$, which is recorded by two separate desktop PC's.

As you can see, there are 3 peaks in the cross-correlation graph. And this is just an example , some recordings gave me about 10 different peaks.


texts: Why do additional harmonics arise after the filter bank?

texts:
I'm trying to implement direct implementation of MDCT filter bank.
So, some additional artefacts have appeared after all of the routines.

The spectrum of the initial signal looks like:

import numpy as np
import matplotlib.pyplot as plt

N = 8 #subbands
L = 2*N

h = np.zeros((L,1)) #skeleton for window function
for n in range(2*N-1):
    h[n] = np.sin((np.pi / (2 * N)) * (n + 0.5)) #window function

w = [30, 100, 300]
a = [1.1, .5, .1]
x = 0

t = np.array([i for i in range(1,3001)])/1000
for idx in range(len(w)):
    x = x + a[idx]*np.sin(2*np.pi*w[idx]*t)

FFT = np.fft.fft(x)
amps = np.abs(FFT) / (len(FFT) / 2)
fs = 1 / (t[1]-t[0])
f = fs*np.array([i for i in range(int(len(x)))]) / len(x)

plt.subplots(1, 1, figsize=(6, 4), dpi=150)
plt.stem(f[:int(len(f)/2)], amps[:int(len(f)/2)])
plt.ylabel('Magnitude of the FFT')
plt.xlabel('Frequencies (Hz)')
plt.grid(True)
plt.show()

# Analysis
H = np.zeros((N,L)) #skeleton for analisis H matrix
pr = (L+len(x)-1) #lengh (number of coloumn) of signal after convolution
Analysis_Mat = np.zeros((N,pr)) #skeleton for matrix after convolution
for k in range(N): #rows
    for n in range(L): #coloumns
        H[k,n] = h[n]*np.cos((np.pi/N)*(k+0.5)*(n+0.5-(N/2))) #analysis H matrix
    Analysis_Mat[k,:] = np.convolve(x,H[k,:]) #convolution

# Downsampling
M = int(Analysis_Mat.shape[1] / N) #number of samles that sould be stay after downsampling
cutmat = int(np.floor(M)*N) #number of rows that sould follow to downsampln block (cut rows that are not fold to N)

Analysis_Mat_DS = np.zeros((N, M))
Analysis_Mat = Analysis_Mat[:, :cutmat]

for k in range(N):
    Analysis_Mat_DS[k, :] = Analysis_Mat[k,::N]

# Upsampling

Analysis_Mat_US = np.zeros((N, Analysis_Mat_DS.shape[1]*N)) #skeleton for signal that should be after upsampling
for n in range(Analysis_Mat_DS.shape[1]):
    Analysis_Mat_US[:,0+N*n] = Analysis_Mat_DS[:, n]

# Synthesis

R = L + Analysis_Mat_US.shape[1] - 1 #legth of signal (number o rows) that should be after convolution
Syntesis_Mat = np.zeros((N,R)) #skeleton for signal that should be after convolution
G = (np.fliplr(H))/(N/2) #Synthesis matrix
for k in range(N):
    Syntesis_Mat[k,:] = np.convolve(Analysis_Mat_US[k,:], G[k,:]) #convolution
y = np.sum(Syntesis_Mat, axis=0)

FFT = np.fft.fft(y)
amps = np.abs(FFT) / (len(FFT) / 2)
fs = 1 / (t[1]-t[0])
f = fs*np.array([i for i in range(int(len(x)))]) / len(x)

plt.subplots(1, 1, figsize=(6, 4), dpi=150)
plt.stem(f[:int(len(f)/2)], amps[:int(len(f)/2)])
plt.ylabel('Magnitude of the FFT')
plt.xlabel('Frequencies (Hz)')
plt.grid(True)
plt.show()

However, the reconstructed spectrum looks like:

Is it OK?
Is there a mistake in my script or I don't understand the theory (or both)?

Can anybody help me with this issue?


texts:
I'm trying to plot some sine waves (code example in plain js here).
When freq is "low" (freq = 10 hz in that case), the plot is quite nice:

The problem is when I increase the freq (try to set var freq = 50 for example):

lots of ripples, it becomes distorted and not so good as plot. If I increment it more, even worse (var freq = 8030 for example is terrible).
When I see those kind of graph on pro systems, they are displayed just fine.
How would you improve it? FFT, splines, whatever? Which is the right approch?
I don't really need accurancy (i.e. for waveform analysis or whatever), just plot it nicely (as in Desmos https://www.desmos.com/calculator/eodkjlywjh, for example).


texts: Algorithm for 1d spline interpolation suitable for 8 bit microcontroler

texts:
What is a concise, fast, down to earth algorithm for doing (or closely approximating) spline interpolation on a 1d continuous stream of data?
(Edit1: The paragraph below equates to saying "the data is uniform in interval." in an awkward way.)
The data is 1d in that sampling on x is at fixed, regular intervals (likely a power of 2 constant) known well ahead of time. This makes y the only variant, which should allow for quite a bit of simplification and precomputation (LUT?).
Here is a graphical example of approximately what we're trying to do. It's our algo applied to a coarsely sampled sine function.
(Edit2: Note that this is only an example of what our algo should do with coarsely sampled sine data, however, the data we would like to process will be random in the set [0,255], most likely 3 or 4 points)

(source: mathworks.com)
Assume high accuracy is not required, but that we must compute the number of results (red dots) between samples (blue circles) in less than 8ms (~120Hz). Also assume the signal processor available is limited in power and instruction set (PIC/AVR), so has only the following relevant instructions and limitations.

(signed + carry) Addition/subtraction instructions.
(unsigned 8x8-to-16) Multiplication instruction. (PIC18, megaAVR only)
Byte wide Boolean instructions (XOR, OR, NOR, AND, NAND, etc.)
Single bit left and right logical shifts. (no barrel shifter)
Can only execute at 2~4 MIPS

Additional Notes:

Responses would preferentially be in pseudocode, so they are more generally applicable.
Cheating is totally OK; it doesn't need to be perfect, just better than liner interpolation.
Bonus points for alternatives that don't need multiplication.
More bonus points for alternatives likely to complete in less than 1ms!

This is for an RGB mood lamp x-mas present for my sister and mommy :3, which I would do myself, but the maths for this are apparently beyond me.

Edit 12-21-2016: Better list formatting


texts: I am new to Kalman and trying to use it for a particle tracking technique for velocity measurments in Multiphase systems. I have particle locations but there are error associated with these locations. Can use of Kalman Filter be useful in reducing error in this case?


texts: Difference between repeated poles and distinct poles?

texts:
An important concept in dsp is marginal stability where we often see the term" repeated roots " or "repeated poles "?
What are they?

Does the term repeated means that two or more poles occur at exactly same coordinates?

Please kindly give answer along with examples for both s domain and z domain?


texts: Average gradient of pixels and sum of absolute diference

texts:
Studying the research paper Multidirectional Scratch Detection and Restoration in Digitized Old Images, E. Ardizzone et al., 2011:

I have several questions regarding this section (5. Restoration),

How can I obtain the average gradient on an area in an image?
How can I obtain a vector that minimizes the sum of absolute differences (SAD) in an image?
Which block-matching approach is being suggested here?
Which factor determines the distance of the blocks from the scratch line?

I could understand the equation (6), but I could not understand how I can implement it practically like in MATLAB or in any programming language.


texts: Calculate harmonics using DFT from real points

texts:
I have a real data of 144 points, when I perform a 144-point DFT on this data, I get $X$ with real and complex values. I want to calculate harmonics using these $X$'s.

The $X[0]$ and $X[72]$, added together and divided by 144, would give me the DC component?
And can I just use the next 71 $X$'s and their conjugates (Euler's identity), to calculate the harmonics? I believe it's $(X[i] + X^*[i])/144 $?

Sorry about that, but I have this problem to solve, I'll try and explain it and what I'm supposed to do with it.
Thanks for being patient.

I have a data set of 144 points. I want to express this in terms of a DC component and it's harmonics. So, I was asked to perform a DFT operation on this and get the harmonics that would represent this data in frequency domain. I performed a 144 point DFT on this data and got 144 $X[k]$'s with real and imaginary parts. As, per my understanding these 144 $X[k]$'s ($X[0]$ being the DC) represent the time domain signal in frequency domain, but I'm still being asked for the DC(apparently ($X[0]+X[72])/2$) and the 71 harmonics, which I'm not sure how to go about.

I apologize if this still doesn't makes sense, but it's what I'm supposed to do, probably you can give me some reference where I can get such concept or anything related.


texts:
I have some difficulty understanding the following question.

I have written a code to plot continuous and discrete version of the chirped signal.
import numpy as np
from scipy.signal import chirp, spectrogram
import matplotlib.pyplot as plt
%matplotlib qt
N = 400
t = np.linspace(0, 2, N)
f0  = 0
f1 = 10
t1= 2
Fs = 8000
n = np.arange(N)
alpha = (f1 - f0)/t1
w = np.cos(2*np.pi*f0*t + np.pi*alpha*t**2)
plt.plot(t, w)
w_d = np.cos(2*np.pi*f0*n/Fs + np.pi*alpha*(n**2)/Fs)
plt.plot(n / (N/2),w_d, 'ro')

which gives,

Although N is completely taken randomly. It gave me a smooth curve for a continuous version. I don't know if my plot is correct or not. However, I am stuck now. what does the question mean? How to find the counts crossing the absicca? For zero-crossing, as far as I understood, the count would vary for N? I'd appreciate it if anyone could help me understand the question.


texts: Group delay of non differentiable unwrapped phase

texts: As we know group delay as a function of frequency, computed as differentiation of  unwrapped phase, but what happen  at those frequencies where unwrapped phase is non differentiable?how to calculate delay at that particular frequency?


texts: Advantage of using nonideal filter in purely digital application

texts:
Say I'd like to bandpass filter a signal with corner frequencies 1 kHz and 2 kHz (i.e., passes signals between 1 kHz and 2 kHz and suppresses other frequencies) using software like MATLAB. Intuition tells me there's two ways to do this:

Use filter design tools to create, say, a 10th order Butterworth filter. Use the filter command with the filter coefficients and input signal. This works, but the filter is nonideal (specifically, there will be some passband ripple, and the frequency cutoffs won't be immediate -- a signal at say, 2.1 kHz, will still be partially passed).
Perform an FFT on the input signal. Set any FFT element corresponding to a frequency outside the 1 kHz-2 kHz range to zero. Inverse FFT the input signal. This will essentially be an ideal filter -- anything between 1kHz-2 kHz is passed with unity gain (no passband ripple), and anything outside 1  kHz-2 KHz is exactly zero (immediate cutoff).

Intuition tells me that #2 is better, because it's an ideal filter. I've learned that an ideal filter can't be built using analog components, but if a certain application is purely digital, why not use the ideal filter?

That being said, I've seen lots of code that uses tools like a nonideal 10th order Butterworth filter -- so what am I overlooking? What's the advantage of using a nonideal filter?


texts: How to deduce the shape of autocorrelation graphically?

texts:
I'm asking for your help in understanding if I have a rectangular function, so the autocorrelation for it will be triangular... and I can prove that graphically but this is time-consuming especially if I have another complicated signal...

So is there any way to know how to deduce quickly the graphic shape of the autocorrelation of any signal?


texts: Why is bandwidth always limited in a real (physical) channel?

texts:
I'm talking about a continuous analog channel.

Why can't it support infinite bandwidth?
Is there a physics reasons for it say for electrical signals?


texts: In the real world, Is there an example of an infinite signal?


texts:
I have signals with a very short rate of ~100e-6 seconds, sampled at 10MHz,so about 1k samples.
I am wondering, is there a practical limit to the usefulness of time-frequency methods concerning a minimal duration of the signal or a certain sampling frequency for that matter?
i am talking about chirp radar signals which fires 1000 linear chirps each 100us long which I sample at 10 MHz. so the whole frame takes only about 100ms


texts: Are the components of a Bandpass signal namely "inphase" and "quadrature" real functions in time domain?

texts:
A band pass signal representation goes by the generalization as
$X(t)=XX(t)*e^{j \cdot2\pi \cdot ft}$
where $f$ be the carrier freq. and $XX(t)$ be the complex envelope.
On further decomposition it boils down to:
$X(t)=XXi(t)\cos(2\pi ft)-XXq(t)\sin(2\pi ft)$

Here is it implicit that $XXi(t)$ and $XXq(t)$ are real functions and their spectra be even symmetric?


texts: How to get an interpolation weight from a mathematical definition

texts:
It was recently explained to me that a "Nearest neighbor" kernel for 1D interpolation can be implemented like this using NumPy

def nearest(delta):
    delta = abs(delta)
    if delta <= 0.5:
        return numpy.asarray([0,1])
    else:
return numpy.asarray([1,0])

Whereas the mathematical definition of nearest neighbor is

$h_{nn}(\delta) =
\begin{cases}
1 & \text{if}& -0.5 \le \delta < 0.5 \\
0 & otherwise \\
\end{cases}
$

Similarly linear interpolation, which in NumPy can be expressed as

def linear(delta):
    delta = abs(delta)
return [delta,1-delta]

But the mathematical definition for it goes

$h_{lin}(\delta) =
\begin{cases}
1-|\delta| &\text{if}&  0 \le |\delta| < 1 \\
0 & \text{if}& 1 \le |\delta|
\end{cases}
$

My question is how to form a kernel with weights out of these mathematical definitions, as the code belonging to them does not paint the same picture that the mathematical definitions paint.


texts: De-Doppler a signal: Audio Timestretch Factor controlled by function f(x)?

texts:
Timestretching usually is constant. Few programs offer bézier handles and such to control the lead-in and lead-out. But that's all. I need better.

So the audio I have was in fact recorded in nonlinearly timestretched form, because the sound source was moving toward the microphone -- geometry, speed, and distances are all known at this point. Basically I want to "de-Doppler" the signal!

The means to take out the Doppler effect is sqrt(1+x^2), x being the ratio of linear to offset distance.

I would now like to apply that function to control the amount of stretching (NOT the stretching algorithm itself!) of the signal at any given point. That function would then serve as kind of an automation curve, only that its done algorithmically instead of by "Touch". The stretch factor could depend on the ratio "actual moment/total length".

Is this at all possible?


texts: Expectation and Auto-correlation of an Output from a Matched Filter.

texts:
[see the attached image]

I recognize the filter in the problem is a matched filter. But I don't understand how this filter actually work, and I am not quite sure about how to calculate the expectations and the auto-correlation in this case. Any help is greatly appropriated.