texts: Get difference equation of a filter from MATLAB transfer function coefficients

texts:
I need to get the difference equation of a specific elliptic filter.
I calculated the transfer function coefficients in MATLAB with:

  %% Low pass design
  n = 10;
  passband_ripple = 1;
  stopband_attenuation = 60;
  f = 4750;
  fs = 44100;
  [B,A] = ellip(n,passband_ripple,stopband_attenuation,f/(0.5*fs),'low');

However, when I use the difference equation, I calculated from the transfer function, the samples go to +-inf. Here is the code where I create the noise input sound and filter it with difference eq.

  xin = rand(1,fs/20);
  xin = xin - mean(xin);
  xin = [xin zeros(1,fs*2)];

  x = xin;
  y = zeros(1, length(x));

  %% Low pass
  b1 = 0.0027;
  b2 = -0.0134;
  b3 = 0.0362;
  b4 = -0.0651;
  b5 = 0.089;
  b6 = -0.0978;

  a1 = -7.953;
  a2 = 29.7485;
  a3 = -68.5714;
  a4 = 107.5494;
  a5 = -119.7162;
  a6 = 95.6848;
  a7 = -54.2043;
  a8 = 20.8340;
  a9 = -4.9109;
  a10 = 0.5401;

  for n = 11:length(y)
      y(n) = b1*x(n) - b2*x(n-1) + b3*x(n-2) - b4*x(n-3) + b5*x(n-4) - ...
      b6*x(n-5) + b5*x(n-6) - b4*x(n-7) + b3*x(n-8)- b2*x(n-9) + ...
      b1*x(n-10) - ...
      a1*y(n-1) - a2*y(n-2) - a3*y(n-3) - a4*y(n-4) - a5*y(n-5) - ...
      a6*y(n-6) - a7*y(n-7) - a8*y(n-8) - a9*y(n-9) - a10*y(n-10);
  end

The filter function does what I want, but I need to be able to filter the sound sample-by-sample.

y = filter(B,A,x);

How can I obtain the right difference equation of my filter?


texts: Matlab: Separate data samples from signal

texts:
I have data samples for an EMG signal of three arm reps. In the data I have 1 signal that is flat then fluctuates in voltage then becomes flat again shows three pulses for each rep. I want to separate the three reps from the signal (chop the signal into 3 parts). I can manually do this by plotting the signal then zooming in and finding the index before the variance in voltage fluctuates and after it flat lines again. I wanted to know if there is an automated way to do this. Some sort of greedy algorithm or a built in Matlab function. Any ideas?

One example:

[12,-23,1,-2,-40,-45,80,85,14,-28,10,-22,-100,-105,80,85,14,-2,
 16,-2,-200,-45,305,85,1,-2,1,-2]

So I would want:

1 [-40,-45,80,85]
2 [-100,-105,80,85]
3 [-200,-45,305,85]

Usually the flat line values vary from zero.


texts:
Lets say I have two objects very close to each other. There is a colour difference between them, but the difference is small. Therefore when finding contours, both the objects are treated like on object. Is there a way to make the difference in colour more significant? As in exponentially increase the small difference in colour? Or if not possible how would I make these objects stand out from each other?

I am using OpenCV please suggest solutions possible in that library.

Thank You


texts: What is the right way to interpolate a 2D grid?

texts:
Let's say I have a 2D grid of temperature measurements in some area and I want to estimate the temperature at some point between the samples. Or at every point, which would basically amount to upscaling an image. Which interpolation method would give the best temperature estimates?

I'm thinking that the "correct" approach would be something like this: assume a prior distribution of temperature in the area, use the measurements to update it with the Bayesian formula and then calculate either expected value or the most likely value at every point.


texts: Understanding the definition of signal's bandwidth

texts:
In class, we have defined bandwidth of a signal as the difference between the maximal and minimal frequency. Also, our tutor explained that the bandwidth indicates how fast the signal is changing.

We have been given an example of a single sinusoid which has zero bandwidth (since it's both the maximum and minimum frequency)

What isn't clear to me is why a rectangular shaped signal has infinite bandwidth.

How can it be explained rigorously and intuitively?


texts: Failed to Get Silence in Audio Streaming

texts:
I want to detect silence in audio streaming. I've been following lot of answers from any website and i feel like know how to do that. but i doubt i'm on the right path coz the result doesn't seems right.

The stream is in mp3 format 16 bit 22050 hz, the data stream divided into small packets. So every packet received I'll check the data if its a header or not to build complete chunk of mp3.
Every complete packets(mp3 chunk) will be decoded into PCM. the PCM data consist of signed value between 0-255. i will split every PCM data into 1.16 seconds of audio to detect silence in that duration. I've read some answer to scale the signals to make it between -1 and 1 for further processing.

My first try is calculating RMS. the PCM data will be splited into 20 miliseconds and calculated to RMS for each number of samples (samples=PCM.length*0.02). Then changes the RMS result to dB with 20*log(rms)/log(10). The result is not as i expected. The more silence in the audio, i get less dB rate (-20dB to -30dB) and get (-50dB to -70dB) for more noise in audio. Then i try to compare it in chart but it doesn't really help at all. I think the result is somehow wrong.

My second try is using FFT. The scaled signals from previous data added with padding to make it 2^n length. The padding used to make signals can be calculated with FFT. After i got the FFT result, then i do A-weighting. I can see some similarity in noise audio. The silence audio gives random chart that can't be compared with noise audio. So i decided to use noise to compare with every a-weighting result.
First i try to compare it with calculating DTW of those signals and make the threshold is not more than 5k in differences, but sometimes i getting audio error(the accuracy is around 90%).
Then i try compare it with another method using correlation ifft(fft(noise) * fft(test)), but i don't know how to use it. I've been comparing same 4 signals and the result doesn't showing any number 1 in that list (the highest is 0.625)

could you tell me what my mistake is? i really confused where my mistake is. actually i want the result like ffmpeg RMS level but i didn't get the same result with ffmpeg in my first try. i can't use ffmpeg library coz my bos ask me to calculate it from memory which is pipe in and pipe out. i've try that and it seems ffmpeg can't be use for that scenario. and tell me if you need some sample of error. I'll provide you that.
thanks


texts: Reconstructing Signal From Its Cyclic Autocorrelation

texts:
Can a signal be reconstructed from its cyclic autocorrelation? Specifically, if we know

$$
R^{\alpha}(\tau) = \int{x(t)x^{\ast}(t-\tau)e^{-j2\pi\alpha t}\mathrm{d}t},
$$

can we reconstruct $x(t)\in\mathbb{C}^{N}$? It seems like the answer should be yes, but it's not obvious to me how to do it. I know there has been work on reconstructing a signal from its autocorrelation, and these algorithms of course require some other knowledge about the signal (e.g., Discrete signal reconstruction from its autocorrelation function and one sample). I guess I am hoping knowing the cyclic autocorrelation will at least relax the required assumptions. Moreover, how does the problem change if we also had the conjugate version:

$$
\bar{R}^{\alpha}(\tau) = \int{x(t)x(t-\tau)e^{-j2\pi\alpha t}\mathrm{d}t}
$$

Edit: Let me try rewording the question. I understand that autocorrelation is not an isomorphism. Is this also true of cyclic autocorrelations? Specifically, can we recover $x\in\mathbb{C}^N$ from

$$
R(\alpha,m)=\sum_{n=0}^{N-1}x(n)x^{\ast}(n-m)e^{-j2\pi n \alpha}.
$$

If we only have knowledge for $\alpha=0$, then we have just the regular autocorrelation and the answer is in general no. However, does knowing $R$ for $\alpha=0,1/N,2/N,\dots$ enable one to recover $x$ or at least reduce the required assumptions necessary for doing so? If not, what if the conjugate version is also available? Can we recover $x$ then?


texts:
I'm trying to filter some motion noise from an ECG signal. To do that I'll try to implement an adaptive filter using the LMS algorithm.

I think I have to calculate the MSE of this:

E[e^2 ] = E[(s + n )^2 ]+ 2E[(s + n)X ]W^T + WE[X^T X ]W^T

in which s+n is the noised ECG,

X (n) = [ 1, Accx (n), Accy (n), Accz (n)] are the values from the accelerometer,

W = [w0, w1, w2, w3] the coefficients.

Wt +1 = Wt + 2αeX,

I'm trying to understand how to calculate alfa and minimise the MSE.

http://www.mathworks.com/help/dsp/ref/adaptfilt.lms.html

can anyone give me a simpler example? I'm kind of lost here.


texts:
I am simulating (in Matlab) the transmit portion of a comm system.  I generate random data, then I want to modulate it to BPSK, pulse shape the data, and then potentially shift it up in frequency by a carrier.  My question is regarding the pulse shaping portion.

Atleast in matlab, to generate the root raised cosine filter using the rcosdesign function, you need to know the roll-off factor  ($\beta$) , as well as the samples-per-symbol ($sps$) and the filter span ($span$).  The samples per symbol will be dictated by the sampling frequency $f_s$ as well as the desired length of the generated signal with a certain datarate $R$.  For example, if I have a BPSK signal with a sampling frequency of 8192 Hz, and I want to generate 1 second of data at a datarate of 256 Hz, then

$$
sps = f_s/R
$$

Now we've established 2 of the inputs to creating the filter are set, the rolloff factor $\beta$ and the samples per symbol, $sps$.  Now, my understanding is that the filter span is kind of arbitrary.  I know that a longer filter span creates a longer group delay, but what are the other tradeoffs to having a longer filter span versus shorter?


texts: Is the Wiener filter operating at the Cramér–Rao bound?

texts:
I have been told (Wikipedia agrees) that the Wiener filter is optimal when signal and (additive) noise are WSS. Optimal in the sense that it minimizes the mean-square error.
The Cramér–Rao bound is the lower bound on the variance of an unbiased estimator of a deterministic parameter.
Does that mean that the Wiener filter operating at the Cramér–Rao bound?


texts: What are Chunks when recording a voice signal

texts:
I am having trouble understanding the example program for audio recording ,
it comes with PyAudio module.

I am using an external Audio Card to record the sound , obviously, the voice should be converted into the digital form to be presented to the computer. But i do not understand how this code converts the input Analog Voice into Digital Form and then saving it.

Here is the code,

import pyaudio
import wave
import sys

CHUNK = 1024
# What is CHUNKS here ?
FORMAT = pyaudio.paInt16
CHANNELS = 2
RATE = 44100
RECORD_SECONDS = 5
WAVE_OUTPUT_FILENAME = "test.wav"

p = pyaudio.PyAudio()

stream = p.open(format=FORMAT,
                channels=CHANNELS,
                rate=RATE,
                input=True,
                frames_per_buffer=CHUNK)
print("start....")

frames = []

for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
    data = stream.read(CHUNK)
    frames.append(data)

print("done...")

stream.stop_stream()
stream.close()
p.terminate()

wf = wave.open(WAVE_OUTPUT_FILENAME, 'wb')
wf.setnchannels(CHANNELS)
wf.setsampwidth(p.get_sample_size(FORMAT))
wf.setframerate(RATE)
wf.writeframes(b''.join(frames))
wf.close()

When the code is executed it record the sounds for 5seconds using an external audio card with a mic.

What I understood from the code is,

The Rate is 44Kbps, which means that 44 kilo bytes are generated from the Analog Audio signals per second.

Then an empty frame is generated frame[] where all the recorded data will be stored.

The next thing is for loop,

for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
    data = stream.read(CHUNK)

the answer to RATE/CHUNK*RECORD_SECONDS is 215.33 which means there will be 215 iterations.

I don't understand why the RATE is divided by CHUNKS and then multiplied by RECORD_SECONDS for the recording of 5 seconds. (RECORD_SECONDS = 5)

Also, What exactly is CHUNK ?


texts: What is the difference between circularly symmetric complex gaussian channel and rayleigh fading?

texts:
What is different between circularly symmetric complex gaussian  channel and rayleigh fading? or are they the same?

I Saw the method of creating circularly symmetric complex gaussian in this website, https://www.researchgate.net/post/How_can_I_generate_circularly_symmetric_complex_gaussian_CSCG_noise ,and it said

  sqrt(var/2) * (randn(1,N)+1i*randn(1,N))

.

On the other hand ,i saw the method of creating rayleigh fading in this matlab website : https://www.mathworks.com/matlabcentral/fileexchange/30477-ber-curve-for-bpsk-in-rayleigh-channel , it said  h =

  1/sqrt(2) * [randn(1,length(tx)) + j*randn(1,length(tx))];

So does it mean that when the variance is 1, circularly symmetric complex gaussian  channel is the same as  rayleigh fading?


texts: What is the algorithm to do a Discrete Hilbert Transform?

texts:
I have a time series data. I have applied Empirical Mode Decomposition and obtained several IMFs.

Now I need to apply a Hilbert Transform to the IMFs to get the instantaneous frequencies.

I have found two papers explaining how to do that:

this one, that says this:

What Dr. Huang proposed for more accurately calculate the instantaneous frequency is to do normalization of IMF. The process is as followed:

Take absolute value of IMF.
Find extrema.
Based on these extrema, construct envelope.
Normalize IMF using the envelope. The FM part of signal becomes almost equal amplitude.
Repeat process 2-4 after the amplitude of normalized IMF retains a straight line with identical value.
Find the instantaneous frequency on the normalized IMF.

and this one that says this:

Notice that the first one says to take the absolute value of IMF and the other one doesn't.

I have tried both but repeating the loop to make all values be contained between 0 and 1 is doing nothing. Every loop the values are the same.

How do I really do that and how do I find the frequency of each IMF?


texts:
I would like to ask a theoretical question concerning the Dirac function. The Fourier Transform of the Dirac function is the value 1 (DC) for every frequency. If we consider the Sampling Theorem, we have to find a maximum frequency in the signal $ \ f_{max} $, so that we can sample with $ \ f_s \ge \ 2f_{max}$. But as we can see from its Fourier Transform, the Dirac function contains every frequency, so we cannot find an appropriate $ f_s $ . My question is, from a theoretical point of view, can the Dirac function be sampled?

Edit: Thank you for your helpful answers guys!


texts: How to calculate resolution of DFT with Hamming/Hann window?

texts:
The frequency resolution of a DFT with a rectangular window of size $N$ is given by $f_s/N$. However, when using other window functions like a Hamming or Hanning window the resolution gets worse. How can I calculate the frequency resolution of a DFT when using a Hamming/Hanning window of size $N$?

I found some numbers regarding "effective noise bandwidth (ENBW)" but am not sure what exactly it means. The numbers were 1 for a rectangular window, 1.5 for Hanning and 1.37 for Hamming window. Does that mean the frequency resolution using a Hanning window is $1.5 \cdot f_s/N$?


texts:
I have a 16bit SAR ADC sampling at 6.4kHz. The Anti-Aliasing Filter is a active 4th order Bessel lowpass (2x LTC6081) with -3dB at 200Hz and -80dB at 3.2kHz.

Now I want to decimate these 6.4kHz to 100Hz. The uC is a 8bit RISC (Atmel AVR) running at 20Mhz.

My first idea was multiple digital bessel biquad filters and downsampling by a factor of 2 but this performs horribly. Then I read about half-band and polyphase filters but I don't have any experience in these topics.

So my primary question is: how should I downsample my 6.4kHz signal to 100Hz?


texts: Can playing a wav file theoretically damage my speakers?

texts:
I've been fooling around with some audio transformations and it's useful for me to listen to them when debugging.  Is it possible to create/play a wav file that will damage my speakers?

I'm working on a macbook pro and using the python module scikits.audiolab to play the sounds.


texts: How to automatically detect blurred document images

texts: The question I have is related to image quality evaluation. Suppose you are supposed to perform OCR (Optical Character Recognition) on a set of document images. Some images, however, are very blurring, and for these images blind deconvolution is needed. However, there are good images that do not need deconvolution. Then my question is: could I find a criterion to automatically separate these images into two groups: blurred image and clear images and then treat them differently when performing OCR?


texts: Understanding the math used to generate a point cloud from a depth image

texts:
In reviewing various algorithms for converting a depth image to a point cloud where the camera
has extrinsic parameters

$c_x, c_y$ the center point of a camera
$f_x = f*m_x$, where $f$ is the focal length of the camera
$f_y = f*m_y$

and $m_x,m_y$ are scale factors relating the width\height of the sensor to the number of pixels
i.e. $m_x $ = (width_of_sensor)/(num_pixels_in_horiztonal_direction) (not 100% sure about this) and similarly for $m_y$.
Then a formula for each $[u,v]$ in pixel space (with its corresponding $z=z_{u,v}$) to be converted to a 3d vector is
$x = ((u - c_x)/f_x) * z$
$y = ((v - c_x)/f_x) * z$
$z = z$
I'm trying to understand the geometry of these equations if they are correct. Any insights appreciated.
Edit: I'm not sure what $c_x, c_y$ are. Are they roughly half the number of pixels in the image sensor? So for example for a depth image of resolution $600 \times 600$, $c_x, c_y = 300$, or rather are they a point in world space, in which case what is the relationship between $u,v$ and $c_x, c_y$ ?


texts:
I want to simulate data from the following model:

$\textbf{z}_k=\textbf{H}\textbf{x}_k+\textbf{v}_k$  $\textbf{v}_k \sim N(\textbf{0},\textbf{R})$

$\textbf{H}$ does not change over time
$\textbf{x}$ is a vector of loadings
$\textbf{R}$ is a diagonal of constants

$\textbf{x}_k=\textbf{F}\textbf{x}_{k-1}+(\textbf{I}-\textbf{F}){\mu} + \textbf{w}_k $   $\textbf{w}_k \sim N(\textbf{0},\textbf{Q})$

$\textbf{I}$ is the identity matrix
$\mu$ is the vector of mean values of $\textbf{x}$
$\textbf{F}$ is diagonal with the AR(1) params which do not change over time
$\textbf{Q}$ is diagonal with the innovation processes for $\textbf{x}$

I have the following code in Matlab

nDates=20000; %number of dates

mats=[1 2 3 4 5 6 7 8 9 10 12 15 20 25 30]'; %maturities
nY=length(mats); %#number of yields

z=zeros(nY,nDates);  %declare vector for yields
x=zeros(3,nDates); %declare vector for factors

R=0.00001;   %standard deviation
I=eye(3); %3*3 identity matrix

v=normrnd(0,R,nY,nDates); %generate residuals
F=[0.9963 0 0; 0 0.9478 0; 0 0 0.774]; %AR(1) matrix
mu=[0.0501; -0.0251;-.0116]; %mean of X
lambda=0.5536;

q = [0.0026^0.5 0 0;0 0.0027^0.5 0; 0 0 0.0035^0.5];
Q=q*q';
rng('default');  % For reproducibility
r = randn(nDates,3);
w= (r*Q)';

B= [ones(nY,1),((1-exp(-lambda*mats))./(lambda*mats)),((1-exp(-lambda*mats))./(lambda*mats))-exp(-lambda*mats)];

X(:,1)=mu;

for t=2:nDates
x(:,t)=F*(x(:,t-1))+(I-F)*mu+w(:,t);
z(:,t)=B*x(:,t)+v(:,t);
end

z(:,1)=[];

It all seems straightforward enough but what test can I so to ensure that it has been implemented correctly?

Ones I have thought of:
Check the correlation of the factors x on their 1 period lags match the values given in matrix F
Check that the variances of v and w are correct
Check that the mean of the simulated variables are correct

I would like to check that the empirical variance of the parameters matches their theoretical equivalent, but I don't know what the theoretical equivalent should be?

Please feel free to suggest further tests that will allow me to know for sure if the implementation is correct.


texts: How can I calculate the Expectation for a particular known vector?

texts:
Let us say I have a row vector $X = [x_1, x_2, x_3, \ldots, x_n]$ and another row vector $Y = [y_1, y_2, y_3, \ldots, y_n]$. I want to check whether the two vectors are statistically independent or not.

Now two vectors are said to be statistically independent if $E(X Y) = E(X)  E(Y)$, where $E$ is the expectation operator.

So how would I calculate the expectation value for each of the concerned vectors?

Thank you for your suggestions - I'm really new to statistics.


texts:
I am using rtl-sdr at an adjustable sample rate to get stuff from a WBFM radio station.  I then filter and demodulate it.  When I sample at 2.4Msps, the sound is 'choppy'.  I tied the sample rate to a slider on gnuradio, so I can vary the software radio sample rate by moving around the slider.  It turns out the pitch of the person talking/singing's voice gets lower as I move the sample rate down, and higher as I move the sample rate up.

Why?  2.4 million samples / second is much greater than 2*20,000Hz.  Why would there be a difference between sampling at 2.4 and 2.8 Msps?  All the frequency content is already covered in both cases...

Also: in gnuradio I can adjust the width of the low pass filter that is placed before the FM demodulator.  I thought that making the filter very narrow would remove the high frequency content of the audio.  Well, it doesn't.  It just makes the sound 'staticy'.  Why?


texts: Designing butterworth filter with pole placement

texts:
I am trying to design a Nth order continous time Butterworth filter in state space, with pole placement technique.

n = varargin{1}
Fc = this.fc
wc = 2*pi*Fc;
T = 1/(2*pi*Fc);
%----------------Finding Poles--------------------------
 p = exp(1i*(pi*(1:2:n-1)/(2*n) + pi/2)); %n is the order of the filter
                 p = [p; conj(p)];
                 p = p(:);
                 if rem(n,2)==1   % n is odd
                     p = [p; -1];
                 end
%----------------Finding Filter co efficents---------------
                 if length(p)<= 2
                    a = conv([1,-p(1)],[1,-p(2)]);
                end
                if length(p) > 2
                    a = conv([1,-p(1)],[1,-p(2)]);
                    for i = 3: length(p)
                        a = conv(a,[1,-p(i)]);
                    end
                end
                a = real(a)
                a = fliplr(a);
%----------------- calculating the matricies%-----------------

            %--------------A matrix---------------
            N = length(a)-1;% order of the matrix
            last_row = (a(1:end-1))*(-1/a(end));% creating the last row

             last_row = (wc.^(N:-1:1)).*last_row;% creating the last row

            k = ones(1,N-1);
            this.x = diag(k,1) ; % diag(vector, k) produce a matrix filled with zero's and k'th off diagonal as vector, in this  case, 1's.
            this.x(end,:) = last_row;% adding the last_row

            %-------------desiging B matrix--------------
            this.y = zeros(N,1); % creates a coloumn matrix with zeros
            this.y(N,:) = 1/(a(end)*wc^n); % writes last elements as a

            %--------------- writing C matrix--------------
            this.z = zeros(1,N); %creates a row matrix
            this.z(1,1) = 1;

            %------------ writing the D matrix----------------

            this.u = [0];

I know the code is bit lengthy, but i have tried everything, but can't figure out where i am going wrong. Because I am not getting any error. The resultant output is just 0. The code is pretty simple, I find the poles, then the filter coefficients. with the filter coefficients. i design A,B;C,D Matrices.
This is a continuous butterworth filter, Can anyone point the mistake i am doing?

edit1:
A low pass PT2 filter:

A low pass 2nd order butterworth filter:


texts: Why does circular symmetry of complex envelopes mean that $\mathbb{E} \Big[\mathbf{x}[i] \mathbf{x}^T[i] \Big] = \mathbf{0}$?

texts:
I am making my way through this paper here, and I understand everything up to equation 16, and I am stuck here.

Long story short, the vector $\mathbf{x}[i]$ here is an $N$ by 1 vector, consisting of the complex envelopes of signal + noise + interference across $N$ different sensors, taken at some time index $i$. We have a desired signal, an interference signal, and AWGN noise at each sensor. Naturally, the desired and interference signals at each sensor are delayed versions of each other, but the AWGN at each sensor is independent of the AWGN at another sensor.

What I do not understand, is the author's claim, where he states that since the desired and/or interference signals are just delayed copies of each other, "the circular symmetry of the complex envelopes thus implies that:"

$$
\mathbb{E} \Big[\mathbf{x}[i] \mathbf{x}^T[i]  \Big] = \mathbf{0}
$$

Circular symmetry of the complex envelopes? What does that mean, and why do we get the above equation from that?


texts: Problems computing the DFT of finite length sequence

texts:
I am having trouble finding the same answer as the solution manual for this sequence.
The problem asks to compute the DFT of
$$
x[n] = \begin{cases}
            1 & \text{for even } n \in \{0\ldots N-1\}
         \\ 0 & \text{for odd } n \in \{0\ldots N-1\}
       \end{cases}
$$
Setting up the DFT equation I agree with the solution manual
$$
  X[k]=\sum_0^{N-1}x[n]W^{kn}_{N}=\sum_0^{\frac{N}{2}-1}W^{k2n}_{N}=\sum_0^{\frac{N}{2}-1}e^{\frac{-j2\pi k(2n)}{N}}
$$
However when apply the closed form geometric series formula $\sum_0^Na^k=\frac{1-a^{N+1}}{1-a}$ I get
$$
\frac{1-e^{-j2\pi k}}{1-e^{\frac{-j2\pi k}{N}}}
$$
while the solution manual says its
$$\frac{1-e^{-j2\pi k}}{1-e^{\frac{-j\pi k}{N}}}$$
Note their lower exponent does not have a 2. Could someone explain where and how the 2 goes away?  Or does someone agree that the provided answer is incorrect?
This is question 8.5c from Oppenheim's Discrete-Time Signal Processing 3e.


texts: Improve performance of speech recognition

texts:
Could anyone recommend good  speech pre-processing/filtering  algorithms/methods that would improve the performance of my speech recognition. I have been primarily been doing image processing related stuff and hence am not much aware of  speech pre-processing methods.

I am currently using Google's Web Speech API for speech to text recognition. This is the current flow of things.

Take input from microphone
Segment the audio using Voice Activity Detection. This is done since Google's API has a limit on the length of the segment it can process.
Send the segment to Google's web speech API (on a thread)
Print the response of the API

Using this I could only get up to 60-65% conversion accuracy. Is there any way to improve this? Would doing a pre-processing on the audio (noise-removal, filtering, etc.) help?


texts: Why is there a dip in this Chebyshev filter response?

texts:
I want to build and simulate a 6 GHz, 4th order lumped-element Chebyshev bandpass filter with a maximum passband ripple of 0.1 dB and a bandwidth of 1 GHz, and as I understand it, there should be 3 lobes in the S11 response, but instead my circuit has a huge drop at the center frequency. I'm not certain what could be causing this or if this is even common, but I'm assuming it has to do with my implementation. The equations I used to get the component values from the normalized lowpass filter coefficients are as follows and the schematic from MWO is also attached;

The elements in the circuit are (from left to right); [C1 = 3.53078138643953e-06 uF, L1 = 0.199281477338676 nH, L2 = 10.3941804502207 nH, C2 = 6.76935843301265e-08 uF, C3 = 5.63670104271011e-06 uF, L3 = 0.124828215212789 nH, L4 = 6.51082585059285 nH, C4 = 1.08069137002873e-07 uF].If anyone can see where I'm going wrong, and would be kind enough to point it out to me that would be greatly appreciated.


texts: What is the purpose of dividing a signal into subbands when compressing an audio file?

texts:
I must implement my own very simple MUSICAM algorithm (basic filtering + Huffman coding) to compress an audio file. But I'm having trouble understanding the interest of dividing a signal into 32 frequency bands before filtering the frequencies. I've read stuff about Nyquist rule but it's not clear for me since I don't have a huge background in signal processing.

Could anyone explain why not just filtering the raw data without dividing into subbands ?

Here is the main reference I use in understanding the concept of MUSICAM.


texts: Python: How to generate log-frequency spectrogram from an audio?

texts:
I want to make a log-frequency spectrogram out of this audio. Later, I need this spectrogram for pitch sequence analysis.
This is a sample sequence I want to achieve:
[ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  0  0  0  0 18 18 19 19 19 19 19 19 19 19 18  0  0  0  0 19 18 18 18 19
 19 19 19 19 19 19 20 20  0  0  0 25 26 26 26 26 26 26 26 26 25  0  0  0
  0  0 26 25 25 25 25 25 25 25 25 25 25 26 26 27 27 27 27 28 28 28 28 28
 28 28 28 28 27 26 27  0 28 28 28 28 28 28 28 28 28 27 26  0  0 26 26 26
 26 26 26 26 26 26 26 26 26 26 26 26 25 25  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0 24 24 24 24 24 24 24 24 24 24 24 24 24 24  0  0  0 25 24
 24 23 23 23 23 24 24 24 24 23 23 23 23 23 22 22 22 23 23 23 23 23 23 23
 22 23 23  0  0  0 23 23 23 22 23 23 23 23 23 23 22 22 23 22 23 22 21 21
 21 21 21 21 21 21 21 21 22  0  0  0 21 21 21 21 21 21 21 21 21 21 21 20
  0  0  0  0 19 19 19 19 19 19]

From what I am able to do now, I can only plot the audio using a linear scale spectrogram, with this code:
def make_spectrogram_b(songname, titles, filename):
  x, sr = librosa.load(songname, sr=None) #sr=None, buat dapet aslinya
  duration = librosa.get_duration(filename=songname)

  print("Audio shape: ", x.shape)
  print("Sample rate: ", sr)
  print("Duration of audio: ", duration)

  # compute stft
  window = np.hanning(window_size) # window size = 1024; hop_length = 256
  stft= librosa.core.spectrum.stft(x, n_fft = window_size, hop_length = hop_length, window = window)
  out = 2 * np.abs(stft) / np.sum(window)

  # plot result
  plt.figure(figsize=(12, 4))
  ax = plt.axes()
  ax.set_axis_off()
  librosa.display.specshow(librosa.amplitude_to_db(out, ref=np.max), y_axis='log', x_axis='time',sr=sr)
  plt.savefig(f'spectrogram_data_B/{titles}/{filename[:-3].replace(".", "")}.png', bbox_inches='tight', transparent=True, pad_inches=0.0)
  plt.clf()

This here is an output sample from code above.
 My testing result is not so satisfying, as it detects too many zero pitch values, and I think I want to change the spectrogram type.
I read from a book source (Muller, Fundamentals of Music Processing, 2015), that if we want to make a spectrogram for music analysis, we must make a log-frequency spectrogram, as quoted:

To emphasize musical or tonal relationships, the frequency axis is often plotted in a logarithmic
fashion, which yields a log-frequency representation. A logarithmic frequency axis also accounts for the fact that human perception of pitch is logarithmic in nature. Finally, in
the case of audio signals, the amplitude values are also often visualized using a logarithmic scale, for example, by using a decibel scale. In this way, small intensity values of perceptual relevance become visible in the image. In the following, if not specified otherwise, we use in our visualizations a linear frequency axis and a logarithmic scale to represent amplitudes. The specific scale is not of importance, but only serves the purpose of enhancing the qualitative properties of the visualization.

In Python, how can I plot this log-frequency spectrogram? Or, is there any better way to 'convert' audio given above to a visual representation for pitch analysis?


texts:
I am a newbie in image processing. in an open source code for convolution neural networks, the author is multiplying a gray scale image from MNIST dataset by 256. Please  can any one tell me why?

ADDITION: the greyscale images have a max value of 0.8681 and a min of -0.1319

Many thanks.


texts:
I am doing my MSc in electrical engineering, and there are not many courses in my university on signal processing. I want to develop my skills in signal processing such that I have a hold in one particular domain.

If anybody could help me and guide me on how to start developing skills it would be helpful.

Skills like which languages to be known. What are the basic things to be familiar with in basic signal processing and image processing? If there is any other domain of signal processing where signal processing engineers are handy, please let me know. Finally what books are helpful in developing the skills in DSP?


texts: Plotting the magnitude response of a filter?

texts:
I have designed a second order butterworth filter.

Sampling frequency: $4\textrm{ kHz}$,   Cut-off frequency: $500\textrm{ Hz}$

$$e = \tan\left(\frac{\pi \times f_c}{f_s}\right) = \tan\left(\frac{500\pi}{4000}\right) = 0.41421$$

Difference equation:

$$y[n] = 0.0976\cdot x[n] + 0.1952\cdot x[n-1] + 0.0976\cdot x[n-2] + 0.9429\cdot y[n-1] - 0.3334\cdot y[n-2]$$

How can I plot the magnitude response of this filter to prove that cut-off occurs at $500\textrm{ Hz}$?


texts: Transposed FIR filter structure to compensate for extra delays on addition

texts:
The problem I'm having concerns an add operation in verilog which is part of a transposed FIR filter.

add(a,delay_line[i],delay_line[i+1]);

a is an input from multiplying a cofficient with the current input sample.
delay_line[i] is an input from the delay line at index i
delay_line[i+1] is an output to the delay line at index i+1

The problem I have is that the add operation takes $3$ clock cycles. My new a input arrives every $1$ clock cycle, so as you can see ideally
I would need the add operation to take one clock cycle, then the
inputs to the add are at the same rate, and can update the delay_line
correctly. I'm sure doing an add in one clock cycle is fine in practice,
but in the case where it is not, then how can this be resolved ?.

My first approach was to draw out the filter structure and add more
delay taps on the ouputs of adds. The normal FIR transpose filter has
single $z^{-1}$ delays after each addition, so I placed some $z^{-3}$ delays to compensate for the add delays. doing that wrecks the frequency response so its not that simple :)

I initially thought this was a Verilog / FPGA problem, but I think this
is more of a DSP algorithm problem now so I posted it to this forum. If someone can point me to some literature on this that would be
appreciated.

Here's my attempt to model this in MATLAB (sorry if you don't have this tool)

%this is an fir transpose filter
fre = ( [0 1500 1700 4000]/4000) ;
msk = [1 1 0 0]; %mask
b = firpm(3,fre,msk); %simple 4 coefficient low pass

%convert filter to transposed form
Hd = dfilt.dffirt(b);
b = cell2mat(Hd.coefficients);

bd2 = [b(1) 0 0 b(2) 0 0 b(3) 0 0 b(4)]; %never this simple!

fvtool(b,1,bd2,1);%compare frequency response


texts:
This is what I know about the basic sampling process

I am told that fbar(t) (i.e the sampled signal) is defined for all t. I don't understand that. If n is a finite integer (..,-3, -2, -1, 0, 1...), then clearly it won't be defined n=1.5. So how is it defined for all t? Also I thought that sampled signals are always discrete and I understand that continious signals are defined for all t, it makes no sense that this is continious.


texts: Fast Calculation of Legendre moments in images using kernel-separability

texts:
I am implementing zero order approximation method to calculate Legendre moment of an image. One simple way to calculate moments is directly using loops and calculate the double summation as given in the image-1. Then I read everywhere that kernels are separable, so using this we can do a faster implementation. Someone please elaborate how separability of kernels reduces time complexity in this case.

Later in the paper - the author says the same things (a little difference is that - he has an exact calculation, but still kernel separability is valid in both cases) about a faster implementation as given in the image - 2.

I will be very grateful, if someone can explain this?


texts: Frequency response of FM modulation/demodulation chain with phase derivative demodulation

texts:
Frequency modulating a carrier by white noise and then demodulating the complex signal using discrete derivative of phase it appears that the discriminator is acting as a low-pass filter.

How do I calculate the rolloff vs modulating frequency for a given phase-based discriminator's impulse response?

I know about Octave's freqz function, but applying it to the discriminator's impulse response gives the frequency response of the discriminator, which in my understanding would translate to amplitude distortion of the demodulated signal, so this is not directly the frequency response of the modulate->demodulate chain which I'm looking for.

pkg load signal;
orig_mono_240Khz = randn(1,150000);
orig_mono_240Khz = orig_mono_240Khz - mean(orig_mono_240Khz);
max_val = max([abs(max(orig_mono_240Khz)), abs(min(orig_mono_240Khz))])*1.2;

subplot(2,1,1);
plot(abs(fftshift(fft(orig_mono_240Khz))));
axis("tight");

phase_changes = j*75000*2*pi.*orig_mono_240Khz/max_val;
signal = resample(exp(cumsum(phase_changes/240000)), 10, 1);
signal = signal - min(signal);
signal = round(signal/max([max(imag(signal)),max(real(signal))])*255);

signal = signal - mean(signal);
rcv_240 = decimate(signal, 10);
unwrp_phase = unwrap(angle(rcv_240));
unwrp_phase = unwrp_phase - mean(unwrp_phase);
disc_resp = [1, 0, -1];
phase_drv = conv(unwrp_phase, disc_resp)(100:end-100);
phase_drv = phase_drv-mean(phase_drv);

subplot(2,1,2);
plot(abs(fftshift(fft(phase_drv))));
axis("tight");


texts: Creating GNU radio flowgraph for continuous switching between transfer and receive mode

texts:
I am having a hackrf-one device which can either transmit or receive signals. I am having a GNU radio flowchart which has a selector block which allows manual switching.

I am not able to switch it at fast rate manually. I am curious to know that is it possible to control the selector block manually so that the reception and transmission can switch at faster rate so i can receive the signal which i am transmitting.


texts: Fourier Transform & Convolution to get response - confused

texts:
in the past couple of years I've always been using Laplace transform for system responses, but today tried to do it with Fourier transform and failed. What am I doing wrong?

If we have a signal $x(t) = 1 + \cos(\frac{4\pi}{3}t)$ and system impulse response $h(t) = e^{-3t}u(t)$, getting $y(t)$ with $y(t) = \mathscr{F}^{-1}(X(j\omega)\cdot H(j\omega))$ should be easy, right?

Well, if we have:

$\displaystyle H(j\omega) = F(e^{-3t}u(t)) = \frac{i}{\sqrt{2\pi}(\omega + 3i)}$

$\displaystyle X(j\omega) = F(1 + \cos{(\frac{4\pi}{3} t)} = 3\sqrt{\frac{\pi}{2}} \delta(4 \pi-3 \omega)+\sqrt{2 \pi} \delta(\omega)+3 \sqrt{\frac{\pi}{2}} \delta(3 \omega+4 \pi)$

The product is an ugly mess, even if we recognize that only the values at frequencies $0$, $4\pi/3$, and $-4\pi/3$ exist, and once we get the inverse Fourier transform, the time solution is complex, which cannot be right, right? :)

What am I doing wrong?


texts: Finding recorded sub-signal in the original audio signal

texts:
I have the original record of a song, and some fragment of the same song recorded by a microphone. I need to find the time where the fragment is appearing in the original.
I have tried to use cross-correlation, the results are good but sometimes it gives a completely different part of the song that is similar (same notes but not the same words).
How can I solve this problem ? Could the cross-correlation value be higher than the time when the sub-signal is found ?
Different approaches are welcome also.
Here is my code for the correlation implementation.

#compute correlation
correlation = correlate(reference,signal,mode='valid',method='direct')
#get maximum index
index = np.argmax(corr)
#extract from the reference
extract = reference[index:index + len(sec)]

You can see the highest peak in the beginning, but I know it should be later. How can I find the real signal in this case ?

Thanks !


texts: Image zoom out using Gaussian blur and downsampling

texts:
I am trying to implement a simple zoom out algorithm using two steps:

The image is blurred using gaussian convolution of standard deviation $0.8\sqrt{{\tt zoomfactor}^2-1}$ as we saw in class.
The image is then downsampled using a simple downsampling of zoom_factor step
The following code is supposed to do the trick.

In class, we saw that zoom out using only downsampling should give worst results than in the case of zoom out combined with gaussian blur. However for my case, this was the opposite which is really weird.

 close all;
 clear all;
 I = double(imread('put an image path here'))/255;

 zoom_factor=2;
 %I=rgb2gray(I);
 sigma=0.8*sqrt(zoom_factor*zoom_factor-1);
 [M,N,s]=size(I);
 f1=-fix(M/2):ceil(M/2)-1;
 f2=-fix(N/2):ceil(N/2)-1;

 [f2,f1]=meshgrid(f2,f1);
 Y=exp(-2*pi*pi*sigma*sigma*((f1/M).^2+(f2/N).^2));
 If=fft2(I);
 If=fftshift(If);
 for i=1:3
   If(:,:,i)=If(:,:,i).*Y;
 end
 newI=ifftshift(If);
 newI=ifft2(newI);

 newI=newI(1:zoom_factor:end,1:zoom_factor:end,:)

 imshow(newI)


texts: Is the common FIR Filter delay line (see pic) also valid on purely logical FPGA implementations?

texts:
Lets say we want implement simply the mean with constant coefficients. So all b's are set to e.g. b=0.5. In this case I do not understand the delay line in the picture below.

I mean I understand it for von Neumann, Harvard architectures.., but how about if I just calculate everything on parallel FPGA? Shouldn't be possible to implement the whole structure without any delays? (..because it is all parallel and there is no feedback) or is it in practice just not possible?

http://en.wikipedia.org/wiki/Finite_impulse_response#/media/File:FIR_Filter.svg

(What I would understand is that there is maybe a delay at the end when the summation takes place)


texts: How do we calculate Power Spectrum Density (PSD) which is given in dB/Hz and not just dB?

texts:
I found many definitions for this quantity, PSD. In a document on signal analysis by National Instruments say that it is the amount of power in a unit bandwidth. Let's say P(f) is the power of the signal contained in a given frequency 'f', then the PSD at frequency $f_n$ is given by $\frac{\Sigma_{i-3}^{i+3}P(f_n)}{f_{i+3}-f_{i-3}}$ where $n \epsilon Z$ and has the units W/Hz or dB/Hz.

I came to know that Power spectrum gives the power contained in a particular frequency (W vs Hz) whereas PSD gives power contained in a unit bandwidth of the signal (W/Hz vs Hz). I realize that they are different but I am unable to figure out how to calculate it. One more thing to add, Gaussian white noise has its power equally spread over all the frequencies and so we expect it to have a constant value in the power spectrum. But we see that it has spikes at some particular values and it varies very much from the constant we expect while in a source I found online PSD remains constant over all frequencies for Gaussian white noise. This is shown in figure below.

I have also found that PSD is mostly used to study all types of noise like white noise, 1/f noise etc. My interest mainly lies in the difference of units between Powerspectrum and PSD which are given by W and W/Hz respectively. Please help me understand this concept if anybody comes across it.


texts: What is the DFT of DFT of discrete signal

texts: What is the discrete fourier transform of the discrete fourier transform of any discrete time signal. Is result same signal? How?


texts:
I was working with message passing blocks (like es_trigger_edge, es_sink and more).
When I connected (using msg_connect) other block in c++ that I wrote, I have noticed that simple processing of the message passing through causes the CPU to load quite drastically.
But, when i used python block i didn't have this problem (I had other problems :)),
So I think I'm missing something important.


texts: Does the Fast Wavelet Transform produce the same coefficient as the Discrete Wavelet Transform?

texts: Does the Fast Wavelet Transform(FWT) produce the same coefficients as the Discrete Wavelet Transform(DWT) if configured for the same depths? Or is the the FWT just an approximation of the DWT?


texts: Fundamental misunderstanding of windowing

texts:
I believe I have a fundamental misunderstanding of windowing. I thought that windowing is applied to a signal with pointwise multiplication, then some sort of filtering is then done to it (say using convolution in the time domain). However, looking back at my previous assignment, I seem to have applied the window to the filter itself(might be remembering incorrectly), which effectively means I convolved the signal with the window? Also, for images, it seems that I need to apply the window to a 2d kernel and not the image itself. This kernel then gets convolved with the image which means that I am once again convolving the image with the window, instead of pointwise multiplying. I understand that windowing is done to prevent spectral leakage/aliasing, but I dont quite get why it is being applied to the filter than the image itself.

If anyone could correct my misunderstanding I would be very grateful.


texts: Spatial moving average in array antenna processing

texts:
In array processing the theoretical spatial response of the array is calculated by doing DFT on sensor locations for some assumed wave number (1/wavelength). The figures below shows the sensor locations of an array (example is from seismology) and the array responses for different wavelength hitting the array. The mainlobe, the blob in the middle of response plots, says something about the resolution while the sidelobes tells of aliasing effects this array configuration gives for assumed wavelength.

In Fiber optic sensing, although every point in the fiber can be considered for measurement, in reality the strain on fiber (because of incoming wave field)  is measured and averaged (moving average) along some length of the fiber, so called gauge length. In general, the averaging is equivalent to multiplying the array response with a low pass filter response in spatial frequency domain (wavenumber domain). The filter response for example in 2D is like:
$$
    \frac{\sin(2{\pi}{k_x L}) \sin(2{\pi}{k_y L})}{L\sin({\pi}{k_x})L\sin({\pi}{k_y})}
$$

where $L$ is the length of the averaging filter and $k_x, k_y$ are the wavenumber components (spatial frequency components). But the averaging examples above assumes regular sampling in either 1D or 2D, whereas the averaging in fiber optics is done along the path of sensor locations. So, in the sensor plot above, the averaging is 1D while the spatial sampling is obviously 2D (in $k_x$ and $k_y$, i.e. $x$ and $y$ components of wavenumber).
I made some synthetic data and performed beamforming to see how the averaging effect is for different wavelengths. But, I need to do this in spatial domain, i.e. find the theoretical array response and multiply it with the response of a moving average filter, where averaging is done along sensor path.
So, my question is: how to calculate the response of the moving average?
Thank you in advance


texts: What's the similarities and differences between Wigner transform and wavelet transform?

texts:
Wigner transform and continuous wavelet transform are both some kind of time-frequency representation of a signal.
What are the similarities and differences between them? Could you give some comparison between them? Let's restrict ourselves in 1D signal at the moment.


texts: How to filter noise from continuous path in image

texts:
Problem statement

In the dataset above, you can notice a continuous path and then a noise-band at around 320-400 on the x-axis. I want to be able to apply a filter that amplifies the continuous path while suppressing the noise-band. The data inside the noise-band is useless.

Attempt 1

I am new to signal processing, but I took a shot at designing my own filter. I tried designing a minimum filter. I defined it to make the output signal $ S $ the minimum of the next $ w $ terms in the input signal $ F $: $ S[n] = \min \{F[x] | x \in \{n, n+1, n+2, ..., n+w\}\} $

for i in range(len(original_signal) - window_len):
    filtered_signal[i:i+window_len] = min(original_signal[i:i+window_len])

It represents each pixel by the lowest pixel in the area which should make the noise go away because the noise band will have a lot of low values in between a lot of high values. It takes out too much of the signal as well.

This is an image showing what the filter took away (the difference between the original signal and the filtered signal)

Attempt 2

A median filter will not work, because the median of the noise-band is relatively high (1e4):

whereas the median of empty space is not that high (1e3):

So the median filter would not reduce the noise to be even with empty space.

Do you know any filters I could try?

Additional Data

Histogram of desired signal:

All histograms normalized and put on the same axes:


texts: Is it possible to generate Hamming code with minimum distance $d$ which can correct more than 1 bit error?

texts:
In literature, all I am seeing is Hamming code with minimum Hamming distance $3$. Theoretically a Hamming code with minimum distance d can detect $d-1$ errors and can correct $(d-1)/2$ error. So minimum distance of $3$ can detect $2$ errors and correct $1$ error.

So here is my question: Is it possible to generate a code with minimum Hamming distance of (let's say) $5$ which can correct $2$-bit errors in each block of coded bits?

If yes, could anyone please give me suggestion about that?

If not, then please give me an explanation.


texts:
According to the Paley-Wiener criterion, a system is causal if satisfies:

$$\int\limits_{-\infty }^{+\infty }{\frac{\ln (|H(f)|)}{1+{{f}^{2}}}}df<\infty$$

So I want to know

This equation is related to LTI system only?
Only casual system is realizable?for linear and non-linear systems?
Which kinds of systems can we create?(only Paley-Wiener criterion can be enough?)


texts: Pearson correlation coefficient vs mutual information

texts:
I recently computed both the pearson coefficient and the mutual information of two time series as well as the corresponding error bars, and what I find out in the end is a significant level of correlation for the pearson r (0.25, with the upper limit of the 95% confidence interval being 0.114), whereas the mutual information shows no significant correlations.

Is this possible or did I make a mistake ?
I understand pretty well what pearson's r is but I find the concept of mutual information a bit harder to grasp and it can be that I misunderstood something.


texts:
Consider multiple images that share a common feature. More specifically, they are the sum of multiple features, and I know that one of those is present in all. What is the best method to extract this common ground?
My naive approach would be to take the minimum over all images for every pixel, but I wonder whether it can be done better.

I am mostly working in Fourier space, but this being a local problem, I expect it to be solved better in real space? A solution in Fourier space would be ideal, though.

For an illustration, see the following images:

What if Goofy was translated and/or rotated?


texts: Mean of input based on acf of filtered output

texts:
By looking at the ACF of the output of a filter whose impulse response I know, can I approximate the mean of the input?

For example

x = randn(1000,1);
y1=filter(ones(5,1),[1],x);
y2=filter(ones(5,1),[1],x + 0.2);
res1 = xcorr(y1,'unbiased');
res2 = xcorr(y2,'unbiased');

When I plot res1 and res2 it is clear res2 has an identical shape but shifted up. This makes me think that just by looking at the ACF I can tell which process has more 'DC' but is there a way to quantify it?


texts: Formula for Bessel low-pass filter coefficients

texts:
When I am filtering a signal in python, it has a built in function to generate bessel filter coefficients given a cutoff ratio and a filter order (number of poles). I am trying to translate this to C code, but I cannot seem to find the formula that is used to calculate the coefficients.

Can someone point me to the formulas to calculate the IIR filter coefficients given a cutoff ratio and order?

I am trying to read the scipy source on github but I am having a very hard time of it...


texts: How to efficiently calculate only the low coefficients of a zero-padded FFT

texts:
I've got an algorithm that zero pads a sequence to 4N, does an FFT, and only uses the lowest frequency N points out of the generated 4N.

This seems like a lot of wasted work, any ideas how this can be done faster?


texts: Is the basis of the sparse signal assumed known in compressed sensing?

texts:
I'm new to compressed sensing, and am a little confused about the assumption of the basis matrix $\Psi$.  Is $\Psi$ assumed known in compressed sensing?

Specifically, suppose that a signal $x$ is sparse in some basis, say $\Psi$, i.e. $x=\Psi\alpha$, where $\alpha$ is $k$-sparse, i.e. $\|\alpha\|_0\le k$.  My understanding is that in compressed sensing, we store $y=Ax$, where $A$ is the $m\times n$ sensing matrix with RIP.  Later on, we can recover $x$ theoretically by

$$\hat\alpha=\underset{\alpha\::\:y=A\Psi\alpha}{\arg\min}\|\alpha\|_0,$$

(or use $l$-1 norm) and then $\hat x=\Psi\hat\alpha.$

In doing so, we need to know $\Psi$ to recover $x$, don't we?  But if we know $\Psi$, why don't we just measure $z=\Psi^{-1}x=\Psi^{-1}\Psi\alpha=\alpha$, and store the value and index of the non-zero components of $z$ ($=\alpha$)?  Can the dimensions ($m$) of $y$ be less than $2k$?  Or is it that $z$ won't be perfectly equal to $\alpha$ due to noise, and we don't want to do thresholding on the fly?  Or some other reason?

I'm confused about the rationale of compressed sensing, and would appreciate any pointers, comments, and clarifications.  Thanks a lot!


texts:
I'm a french student and I have to make a scientific project in group, so we decided to work on how automatically turn sheets music when a musician reaches the end of the page. For that I'm working on pitch detection but my math/physics level doesn't allow me to understand everything (I'm 17 years old).

Now, I'm doing the FFT of a wave signal (microphone input) and get the peak. This works well with sinusoids and wind instruments but when I try whith a violin the results are weird (because there are a lot of harmonics stronger than my fundamental)...

So I have some questions :

Can you recommend me a simple website explaining basic DSP or algorithms to do pitch detection ?
Do you think that just a FFT will do the job ?
There are a lot of methods to do what i want to achieve (autocorrelation, FFT, Cepstrum...), which is the best in my case ?

I really want to understand what I'm doing, that's why I'm coding everything instead of just use a library, but for the moment I'm really lost...

Thank you! (and excuse me for my poor english ^^)


texts: video resources for learning image processing?(for EE background)

texts:
Please kindly recommend some youtube links for learning image processing. I have EE background and i want to study/learn image processing in this way

Also i have confusion that Computer science professionals also do teach/study image processing but that image processing is somewhat different from the one that EE people teach/study


texts: Find camera pose by knowing triangle coordinates

texts:
I have this concrete problem where I know the 3 coordinates of a triangle $ABC$ in a $xyz$ frame reference and I need to estimate the position $P$ and direction $\vec{n_{p}}$ of a stationary camera with known focal distance $f$. For this purpose I have the image from the camera where I know $H$, $W$ and the respective 2D positions $a$, $b$, $c$ in pixels.

My questions:

Is this possible by knowing only 3 points $ABC$?
Does make difference if triangle $ABC$ is equilateral or not since I know the map $ABC => abc$?
How to estimate $P$ and $\vec{n_{P}}$ relative to frame $xyz$?

I read already Camera Pose Estimation from Vanishing Points, Step by Step Camera Pose Estimation for Visual Tracking and Planar Markers and How to compute camera pose from Homography matrix? but sincerily I couldn't realize how to solve my problem.


texts: Speech recognition algorithm in GNU Octave to convert speech to text (and vice versa)

texts:
Recently, I am building an application in Octave and I saw some algorithms about speech recognition spread around the net, but most of them were implemented in Python, MATLAB, C++ and no one in Octave. I already came to test the MATLAB one in the Octave IDE but there are lots of packages that exists in MATLAB and there aren't in Octave, I have tried already to look for at the net in order to find anyway to run extensions to run Python functions in Octave for example but I've only found simple ones to return list of numbers, these things...
Is there any possibility to implement a speech recognition integrated with google voice that converts speechs in text in GNU Octave?


texts:
I try to filter an image using this equation:
$$
I_f(x_j, y_k, t)=\sum_{p=-q}^q\sum_{r=-q}^q I(x_{j+p}, x_{k+r}, t)
$$

where $I$ is the original image in grayscale, $I_f$ is the filtered image and $d$ is the subset. Here, I chose 3. However, I'm confused what will happen when $j=1$, $k=1$, and $p=-q$, $r=-q$.

Looping in MATLAB will result in error:

I=rand(240, 320);
[m, n]=size(I);
I_f=zeros(m,n);
d=3;
q=(d-1)/2;
for x=q:n
    for y=q:m
        for r=-q:q
            for p=-q:q
                I_f(x,y)=I(x+p,y+r);
            end
        end
    end
end

How to do it correclty in matlab? Is it possible to do it without looping? Thanks in advnace!


texts: channel impulse and frequency response: H?

texts:
The received signal of MIMO system is given as:
$$ y= H \cdot x +n, (*)$$
where $H$ is a channel matrix.

$h(\tau_l)$ is a channel impulse response. It is determined as a matrix $\in \mathbb{C}^{Nr \times Nt}$, $\tau_l$ is delay.

My question: $H$ in (*) is a frequency response, fourier transform of   $h(\tau_l)$, isnt't?


texts: How to calculate the gain in a bivariate fft in R?

texts:
In Statistica gain is defined as follows:

  Gain. The gain value is computed by dividing the cross-amplitude value
  by the spectrum density estimates for one of the two series in the
  analysis. Consequently, two gain values are computed, which can be
  interpreted as the standard least squares regression coefficients for
  the respective frequencies.

However, spec.pgram (the engine behind spectrum) in R does not return the cross-amplitude value (as far as I can tell).  How can I calculate gain for these signals?

Example:

makewave <- function(freq,phase,amp,Nsamples=length*samplerate,samplerate,time=Nsamples/samplerate,as.time.series=TRUE) {
  time <- Nsamples/samplerate
  phase <- phase*(2*pi)/180
  wavetimes <- seq(0+phase,time*freq*pi*2+phase,length.out=Nsamples)
  #plot(1:samples/samplerate,amp*sin(wavetimes),type="l",xlab="Time")
  if (as.time.series) {res <- ts(amp*sin(wavetimes),deltat=1/samplerate)} else {res <- amp*sin(wavetimes)}
  return(res)
}
signal1 <- makewave(15,30,1,180,60)
signal2 <- makewave(15.01,60,1,180,60)
signal.union <- ts.union(signal1,signal2)
sp <- spectrum(signal.union,span=c(3),taper=0)
plot(sp)
plot(sp,plot.type="phase")
plot(sp,plot.type="coh")


texts: Input and Output SNR Calculation in AM receiver input and output

texts:
For the FOM calculation of an AM receiver, the input and output signal powers are respectively,
$$P_{in} = \frac{A_c^2(1 + k_a^2P_m)}{2}$$
$$P_{out} = k_a^2A_c^2P_m$$
These are obtained from the input and output expressions for the receiver as given by,
$$x(t) = A_c(1 + k_am(t))\cos(2\pi f_ct) + n_c(t)\cos(2\pi f_ct) - n_s(t)\sin(2\pi f_ct)$$
$$y(t) = k_aA_cm(t) + n_c(t)$$
where $n_c(t)$ and $n_s(t)$ are the in phase and quadrature components of the complex noise envelope, and the output is obtained via envelope detection. Now, while calculating output noise power, we multiply the PSD of $n_c(t)$ with the bandwidth. So,
$$P_{noise,out} = 2N_0W$$
because $n_c(t)$ has a PSD of $N_0$ spread out in a frequency range of $-W$ to $W$. This gives us the SNR of output as,
$$SNR_{out} = \frac{k_a^2A_c^2P_m}{2N_0W}$$
Now, in the literature, as in Wikipedia for instance, (link below) the input SNR is,
$$SNR_{in} = \frac{A_c^2(1 + k_a^2P_m)}{2N_0W}$$
which corresponds to
$$P_{noise,in} = N_0W$$
Now, it is my understanding that the input noise is a bandpass AWGN noise with a double sided PSD of $\frac{N_0}{2}$ spread out in a frequency range of $f_c - W$ to $f_c + W$ and $-f_c - W$ to $-f_c + W$, which then should correspond to a noise power of $2N_0W$. But it seems that the negative frequencies have not been accounted for while calculating input noise power. Where am I going wrong here?

Link to wikipedia page:
https://en.wikipedia.org/wiki/Signal-to-noise_ratio


texts: Finding Minimum SNR or Minimum Detectable Signal Based on Noise Spectral Density

texts:
Imagine I have an $n$-bit sensor, whose noise spectral density is $V/\sqrt{Hz}$. $V$ is an exemplary measurement unit. The desired measurement bandwidth is also known. Is there a mathematical way to find a minimum detectable signal or minimum SNR or clean (unaffected by noise) sensor bits or ... based on the available noise density?
As far as I understand, minimum detectable signal concept is normally used for sender/receiver systems. However, for my sensor, I would like to know how powerful the measurand quantity has to be, such that it can be distinguished from the noise.
I will very much appreciate it, if you provide references for your answers.


texts: Compute RMS value of non-stationary signal

texts:
I would like to know what is the most accurate way to compute a measured signal's fundamental RMS value as a function of time.

The following is true for the considered signal:

signal represents a voltage
sample rate: fs = 6 kHz
capture time: t = 100 s
fundamental frequency: f1 = 60 Hz
signal is not a pure sinusoidal but has additional frequency content
amplitude of the signals fundamental changes slowly over time


texts:
I read from Wikipedia that signals obey the uncertainty principle. The product of the measured uncertainty in frequency times the measured uncertainty in time is greater than a constant. The localization of both frequency and time is not possible simultaneous.

Can two people play a signal game with signal acquisition? Player A sends a signal to player B. Player B listens to the signal. Player A stops transmitting the signal after a short period of time. Player B must calculate from the limited time signal what the frequency spectra is

Is this a real game in the world?


texts: When does the narrowband model in array processing fail?

texts: I am not sure that I understand well the narrowband assumption in array processing. I wonder how it simplifies the model and what would change if not assumed. Also, what will happen if I apply the model on signals that are not narrowband?


texts:
Im using CC1350 launchpad microprocessor of Texas Instrument that I configure my transmitted packets using SmartRF studio
Im trying to decode my packets for verifying if my transmitted packet is valid or not valid by checking CRC- attaching down my transmitted packets that I configured :

So implicitly Im trying to build a receiver side in software - using matlab.
I have a dongle SDR-RTL that Im sampling my packets with frequency sampling 2048KHZ , the data of the dongle output is a sampled file which it's IQ SAMPLES (RAW DATA) I plot it in my matlab and I see the packets - attaching down a photo of my plot:

Im using Legacy Mode which it's documented here over Texas Instrument Handbook :
enter link description here
as you see in the handbook of texas I have in long range mode (legacy mode) FEC encoder and DSSS blocks so in my case DSSS=8 I must see the packets that Im sampling on the output of my demodulator as terms of 11001100 or 00110011 (in terms of DSSS=8 as attached on texas handbook)
So I built a demodulator (my demodulator is two matched filters using Exponential notation EXP(2PIFC))
on the output of my demodulator I should see the decrypted data of IQ samples (of the Raw Data) but Im not seeing my syncword or preamble data in terms of DSSS=8 (in terms of 11001100 or 00110011 according to table of DSSS ) , could anyone instruct me what should I do? my demodulator I build it by using two matched filters (two exponents) and I used zero crossings for sampling the middle points of zero crossing points, my demodulator is works fine in general Mode (not Long range mode) but in Long range mode my demodulator isn't working fine, any help?
why it's not working fine , because as attached in the photo of smartRF over legacy mode packets I must see on the output of my demodulator the syncword or preamble in terms of DSSS=8 (i.e in terms of 11001100 or 00110011 ).
Could anyone help me or instruct me what should I do in order to succeed to see my syncword or preamble in the output of my demodulator? Im in a doubt if I must see my syncword or preamble of my packets on the output of my demodulator as terms of DSSS=8 (11001100 or 00110011) even before doing decoder or encoder code in matlab ...... must I see in my output of my demodulator the syncword or preamble of my packets in terms of DSS=8(11001100 or 00110011)?
thanks alot
Blockquote


texts: How to Tx CDMA signal physically using an existing point-to-point modem?

texts:
I am working in acoustic communication, I want to test CDMA network capabilities in acoustic environment. I had an acoustic modem that has serial interface, a WHOI micromodem.
How can i serialized coded data from various node in the same channel?


texts: How to correct for "track jumping" for multi-lined time series

texts:
I am creating an algorithm that tries to estimate three time-series that should not cross each other (voice formants). The problem is sometimes the algorithm mistakes one time series with another and the algorithm "jumps" track.  Here the blue line/orange line shares one track and orange line/green line also shares a track.

I was wondering if there is a name for algorithms that would fix these kinds of issue. An algorithm that estimates the best fit lines that minimize the number of sharp jumps for example. If this is a developed field I'd appreciate it if someone can point me to the right direction.


texts: How to Deduce a Linear System's Impulse Response from a Set of Input and Output Signals?

texts:
I want to know how to solve those types of problems.. is it by inspection ?

Consider the linear system below. When the inputs to the system $x_1[n]$, $x_2[n]$ and $x_3[n]$, the responses of the systems are $y_1[n]$, $y_2[n]$ and $y_3[n]$ as shown.

Determine whether the system is time invariant or not. Just your answer.
What is the impulse response?

EDIT:  Assuming a general case where the given inputs don't contain a scaled impulse like $x_2[n]$


texts: Identify the frequency spectrum of an audio file by using matlab

texts:
I need to identify the first few harmonics of a audio file by cancelling the noise of the specific signal.
So I have partitioned the signal into 1024 Hz sample blocks and then took the FFT values of each block separately and finally add all FFT values and plot the waveform. Here is my matlab code. Somebody help me with the code since i cannot add the first block fft into total fft. Get an error. Please somebody help me

clc
clear
[signal,fs] =audioread('Female_55.wav'); %Only for plotting
t=0:1/fs:(length(signal)-1)/fs;          %Only for plotting

fs0=fs/1024                        %number of frequency blocks
t1=length(signal)/fs;              %length of the waveform

s=1;    %starting value
s1=1;   %ending value
c1=1;

wavefft_tot= zeros(1025,1);

for blocks = 1:fs0

[signal1,fs1] =audioread('Female_55.wav',[s*c1 s1*1024]);
wavefft_next=abs(fft(signal1));
wavefft_prv=wavefft_next;

n=length(signal1)-1;  %plotting information
f=0:fs1/n:fs1;        %plotting information

if(s==1) && (c1==1)
    c1=1;
    s=1024;
    s1=s1+1;
    wavefft_tot= wavefft_next;
else
   s1=s1+1;
   c1=c1+1;
   wavefft_tot=wavefft_tot+wavefft_prv;

    end

    end

    figure(1);
    plot(f,wavefft_tot);
    xlabel('Frequency in Hz');
    ylabel('Magnitude');
    title('The Wave FFT');


texts: Can compressed sensing be used instead of intepolation for missing values?

texts:
Consider a signal that is sparse in frequency, but it measured in the time domain, for example (in matlab):

t=linspace(0,4,2^8);
x=sin(2*pi*t)+sin(2*pi*3*sqrt(2)*t+pi/4);

let's say that some fraction of the signal is missing, for example the 120-130 measurement:

x(120:130)= NaN;

using standard interpolation methods (nearest neighbors, spline, etc) will fail to reproduce the missing signal, because these methods do not make use of frequency information of the signal.

I thought that by using compressed sensing (L1-magic for example) and random sampling of the signal outside the missing part I'll be able to recover the signal and the missing part, but looking into it I saw that I still need the full signal to create a measurement vector via y=Ax, where A is the sensing matrix.

If x has NaN values then y will be all NaNs, and the entire process is compromised.

So how can I "compress" my sampling by not having the full signal even though it has the information needed? (assuming the sparsity above)...

I'm aware I may use other techniques but I'm curious how to use CS.


texts: Linear prediction (LPC) of Sine wave samples around maximas

texts:
I think I missed class when this was explained ...

Anyway as part of a bigger project I have to implement a LPC to predict 2-3 future values of a sinusoidal process. I wrote a small Matlab m-file to calculate the predictor coefficients and plot the resulting predicted values.

It does look well so far when moving around $y(x) \cong 0$ but the values around the maxima/minima $y(x) \cong \pm 1$ are completely off (it predicts into the opposite direction):

I was expecting something like this:

I'm using the variables 'phi' to set sort of a operation point and 'unk' to set the starting point from where values should start to get predicted.

Here my m-file:

N = 100;            % samples
fs = 60;            % sample-frequency [Hz]
Ts = 1/fs;          % sample-time [s]
unk = 10;           % future values to predict (mostly used to set x(k) where we start predicting from)
P = 10;             % predictor order
phi = 0;            % default phase-offset
phi = pi/2;

k = 0:1/N:2;
x = sin(2*pi*fs*k*Ts+phi);  % reference values

x_pred = x(1:101);          % intialize with x ...
                            % last 'unk' values will be replaced by the predicted ones
x_kn = x_pred(1:end-unk);   % values we already know (we will try to predict last 'unk' values)

coeffs = aryule(x_kn, P);      % get 'P' filter coeffs using the known input 'x_kn' values

for i = 1:unk
    % x^ = -vec(a) * vec(x_kn)'
    nextValue = -coeffs(2:end) * ...            % P - 1 coefficient vector
                x_pred(length(x_kn)-1+i: ...    % last value we actually know (at pos len(x_kn))
                -1: ...                         % reverse iterate ...
                length(x_kn)-P+i)';             % until P values before len(x_kn)
    x_pred(length(x_kn)+i) = nextValue;         % now we know the value at len(x_kn)+1
end

figure(1);
plot(0:length(x_kn)+unk-1 ,x_pred, 0:length(x)-1, x); grid;
legend('predicted', 'real');


texts: Are there any lossless embedding algorithm to support Haar(DWT) in audio watermarking?

texts: I've been trying to do some research on audio watermarking and what i wanted is to insert string into audio file in the transform domain (DWT), i've read some literature of DWT especially HAAR transform and try to implement it. So far i can embedd and extract watermark from the audio which is wav file beacuse i following this literature but apperently it is not lossless watermarking as i try to implement it. What i want to know is from your logical view is that literature really give a lossless method because it said that they got the analisis which mention it supposed to be lossless (BER equal to zero or near zero) when attacked. Or maybe i'm the one who implement it wrong, my analisis give me result that the proposed method is not lossless(sometime the result is wrong). If there is any embedding method that support haar and which is lossless, because i hardly find the literature that exactly pointing out audio watermarking (that insert some text to audio).If anyone know any lossless embedding method in transform domain please help me?


texts: RoC of Z transform of signal consisting of 3 values

texts:
Let's consider the signal $x[n]=\{x[1],x[2],x[3]\}$

It's $\mathcal{Z}$ transform is $\frac {x[1]}{z}+\frac {x[2]}{z^2}+\frac {x[3]}{z^3}$

My textbook says it converges for all values of $z$ but there is a pole at $z=0$, so which one is true?


texts: Why do we use normalized angular frequency?

texts:
Let's say we want to generate a signal in MATLAB, using sin, representing the A-major chord (440 Hz, 550 Hz, 660 Hz, each with amplitude 0.1).

I know we have to use the normalized angular frequency. But I don't understand why? Could someone kindly expalin this to me, please.


texts: I have to do a presentation about Auto-Tune and its relation to the Fourier transform. What is a good explanation on how does Auto-Tune work?


texts:
There are nice technical definitions in textbooks and wikipedia, but I'm having a hard time understanding what differentiates stationary and non-stationary signals in practice?

Which of the following discrete signals are stationary? why?:

white noise - YES (according to every possible information found)
colored noise - YES (according to
Colored noises: Stationary or non-stationary? )
chirp (sinus with
changing frequency) - ?
sinus - ?
sum of multiple sinuses with different periods and amplitudes - ?
ECG, EEG, PPT and similar - ?
Chaotic system output (mackey-glass, logistic map) - ?
Record of outdoors temperature - ?
Record of forex market currency pair development - ?

Thank you.


texts:
I use scipy.signal.decimate  to downsample a large spectroscopic data-set.
This has worked superbly so far.

yy = decimate(data.row(i).buffer, n2, ftype = "fir")

I used successively to downsample the spectrum to several lower resolutions.

When I moved from scipy 0.17 to 0.18 (actually 0.18.1) during of a python3 port, a problem appeared.
The first downsampling works, on all the rows; but the next one crashes with the following message:

File "..../lib/python3.5/site-packages/scipy/signal/signaltools.py", line 3049, in decimate
raise TypeError("q must be an integer")
TypeError: q must be an integer

I looked to the code, but with no avail. I also tried the new zero_phase=True argument, but did not help. Needless to say that my n2 argument above IS an integer, the second time as well as the first time!

Any pointers ?


texts:
i use  this for complex fft.

Output expected

$fft[3].real= 32$ (peak at 3rd bin)

$fft[61].real= 32$ ((peak at negative frequency pair of 3rd bin))

All other values negligibly small

The input is

$ y.real = \sin (2\pi*3i/64)$ where $i = 0 \to 63$

$ y.imaginary = \sin (0*3i/64)$ where $i = 0 \to 63$ (all zero)

The output i got

$fft[3].imaginary = -32$ (peak at 3rd bin)

$fft[61].imaginary = 32$ ((peak at negative frequency pair of 3rd bin))

All other values negligibly small

This is the first time i am working with complex input fft. Can somebody explain me why i am getting peaks in imaginary part and not in real part of fft. As of my understanding doing Real FFT is nothing but using one half of input as Real part input and other half as Imaginary part input to the Compelx radix 2 FFT algorithm [i remember from John G Proakis textbook]. But i don't understand why this pseudo complex signal produces this kind of output. Also can someone explain me the phase information of this pseudo complex
wave ?


texts:
From what I've read, the spectral properties of the modulation largely depend on the rate of change of the signal phase - "a signal in which the phase changes stepwise has a broad spectrum, well beyond the allocated radio channel". I don't know if I understand this correctly.
An example of the spectrum of QPSK modulation with a random signal causing step changes in the output signal is shown in the figure below:

The book I am using says that "if we assume that the bandwidth of the radio channel allocated to this emission is $B_k = 1/T$, then a significant portion of the signal power (approximately a quarter) is emitted outside the bandwidth of that channel. The appearance of such strong interference in other radio channels is unacceptable (moreover - the band in which 99% of signal power is transmitted is approximately ten times wider than $B_k$)"
What does it result from? Why ~1/4 of the signal power is emitted outside the bandwidth of this channel? How do phase spikes affect this? And why is the bandwidth in which 99% of the signal power is transmitted approximately ten times wider than $B_k$?


texts: Correcting large frequency offsets at very low data rates

texts:
Consider a typical BPSK system running at a baud rate of 1 kbaud/sec - the receiver needs to deal with severe frequency offsets, let's say up to +/- 1 MHz offset and that the frequency offset can change by up to 100 khz/sec.
So we're dealing with a frequency offset that's 1000x the baud_rate. Usually if we were transmitting at a higher baud rate, this would be more like 0.1x the baud rate - and one way to easily handle this frequency offset would be to employ fred harris' band-edge frequency locked loop (fll). However, AFAIK the band-edge fll approach only works very well at 2 samps/symbol. If you try to run it at a much higher # of samples/symbol (e.g. in this example we'd need to run it around ~2000 samples/symbol to have unambiguous view of our maximum frequency offset), the band-edge filters would become extremely narrow and you'd have no hope of exciting them with enough energy and locking onto the signal.
Of course - we could just replace the band-edge filters with much wider ones in hopes of catching the signal and locking onto it - but in doing so we're now allowing in a significant amount of noise. This leads me to another approach - where an FFT loop is used to determine the frequency offset and correct for it. I think of this as a generalization of the idea of swapping out the band-edge fll filters for wider ones. For example, we can consider creating a frequency discriminator/error term by comparing the pos and neg half energy of the FFT. Again, this is plagued by similar noise issues. Additionally, we have to consider complications due to the frequency offset moving (at 100khz/sec) - if we try to use a larger FFT (or average multiple FFTs ala welch/bartlett/etc) we run into the issue of the signal smearing across multiple bins - this makes it difficult to achieve significant processing gain and correct for large frequency offsets in low-snr situations.
So my question is - is there a 'correct' way of dealing with this problem? I've convinced myself using a coarse frequency correction followed by a fine freq correction (say an FFT method followed by a band-edge method) sort-of works, but I'm not really sure this is the best way to go at it (especially given the noise issues inherent in the FFT trade-off I mentioned above).


texts: Get slopes of integer shift of an image in radians

texts:
I tried to reproduce the result from this paper, However I'm stuck at the Section III of the paper which involve integer shift. Two images which their phase correlation $\theta(k_1,k_2)$ was calculated at pixel-accuracy(integer shift). I get their phase component as shown(from -$\pi$ to $\pi$). It can be represented by $\theta(k_1,k_2)=\theta_1(k_1,k_2)-\theta_2(k_1,k_2)$. Since the two images shifted in parallel, it can be approximated as $\theta(k_1,k_2) \approx ak_1+bk_2$  where $a$ and $b$ is the slopes obtained through least square method, $k_1$ and $k_2$ are the value $MxN$ of the image. Plotting $\theta(k_1,k_2)$:

And to get phase component with only sub-pixel accuracy(decimal shift), the phase component is subtracted with integer shift. The equation is $\theta''(k_1,k_2)=\theta(k_1,k_2)-(a'k_1+b'k_2)$

However, as you can see from the second figure, the value of $\theta''(k_1,k_2)$ ranges wider than -$\pi$ to $\pi$(as compared to Figure 4 in the paper), which I suspect the slope($a'$,$b'$) was obtained inaccurately. The paragraph above the Equation (14) mentioned that they used conventional POC to obtain the slopes. However I'm not sure did I understood it correctly.

The code I implemented:

function [ output_args ] = phasecorrlsa( refIm, shifIm )

F=(fft2(double(refIm)));
G=(fft2(double(shifIm)));
[m, n]=size(refIm);
[M,N] = meshgrid(1:m,1:n);
X = [M(:), N(:)];

R=(F.*conj(G))./abs(F.*conj(G));
r=(ifft2(R));
% surf(r)
[apx, bpx, rhat]=lsa(angle(r));
[~,w] = max(r(:));

[del_hat2p, del_hat1p] = ind2sub(size(r),w);
del_hat2p=del_hat2p-1;
del_hat1p=del_hat1p-1;

ap=(del_hat1p*2*pi)/m;  %THE SLOPES,obtained through Eq. (12)
bp=(del_hat2p*2*pi)/n;

theta=angle(F)-angle(G);
R=exp(1j*theta); %E6
theta=atan2(imag(R),real(R)); %E9

[a, b, thetahat]=lsa(theta);
del_hat1_ts=(m/(2*pi))*a;
del_hat2_ts=(n/(2*pi))*b;

thetapp=theta-(ap*M+bp*N);
figure;
surf(theta);
figure;
surf(thetapp);

[app, bpp, thetapphat]=lsa(thetapp);
del_hat1=(m/(2*pi))*app+del_hat1p;
del_hat2=(n/(2*pi))*bpp+del_hat2p;

if del_hat1>n/2, del_hat1=del_hat1p-m; end
if del_hat2>m/2, del_hat2=del_hat1p-n; end

output_args=[del_hat1 del_hat2];

end

function [ a, b, hat ] = lsa( theta )
[m, n]=size(theta);
[M,N] = meshgrid(1:m,1:n);
X = [M(:), N(:)];

B=regress(theta(:), X);

a=B(1);
b=B(2);
hat=reshape(X*B,m,n);

end

Anyone could lighten me what went wrong with my code? Thanks in advance!


texts: 2D Convolution in the Spatial Domain vs Frequency Domain

texts:
Suppose, I have this kernel.
-1, -1, -1,
-1,  9, -1,
-1, -1, -1

Can this kernel be used in a FFT based convolution? How?
What could be the reason of my failure?

Related:

Image Convolution in Frequency Domain.
FFT Convolution - 3x3 Kernel


texts:
I'm new to dsp.stackexchange and digital signal processing. I'm required to develop an app to measure the similarity of two signals.for an example I need to compare breathing of a man whether it is abnormal or normal. I have recorded two breathing signals.

Relative error of two signals

I have compared two signals and obtained correlation of two signal. It gives crazy result (because I'm using Android phone it might contain noise). I would like to calculate the  relative error of two signals. In the above link there is an answer on how to get the relative error. But I can't understand the formula. I have two double arrays containing values after doing FFT. I'm stuck at applying above formula. How can I go through this?


texts: Online boosting (training on the fly) for face detection

texts:
I want to implement an on-the-fly training process for face detection.  I'm not sure but it seems it's called online boosting.

Вefore, I have used OpenCV face detection with my own cascades, but I'm interested in technique which can "learn in process".

For example, If I train a boost classifer with 1000 marked images and then I look at the results on my test set, if I'm satisfied, the algorithm will finish. Otherwise, I will add 100 images more but then not retrain all 1100 images, just add 100 images to existing cascade, is it possible? It's also usefull for training in which can participate many people or computers.

What if I have big in-class difference, for example, frontal faces and rotated faces.  Do I need to train 2 separate classifiers?


texts: Role of Microphone sensitivity in building an SPL meter or dB meter

texts:
In idea of buying an SPL meter i have drawn the specifications out in them very first and important spec is sensitivity of microphone many high priced microphones have their sensitivity arround 50mV/Pa where as some other with less cost have it arround 1V/Pa,
so what actually i want to know is, what is the actual role of Microphone sensitivity ?,for me it looks like the mic sensitivity will help in extending its dynamic range only i think..
for example the low cost meters will have differnt dB ranges

low    30-80dB, Med 50-100dB,High   80-130dB
why so

which microphones will be costly,whether a high sensitive like 12mV/Pa or a low sensitive 40mV/Pa


texts:
In simulating linear frequency modulated (LFM) chirp signals, with the parameters:
Pulse width, bandwidth.

I find the matched filter response to be this:

In particular, there is a increase in the energy over the right side of the plot.

My question is, is this normal, or should the tapering be going down?


texts: How is the matrix $\mathbf R_x$ not Toeplitz in case of a signal missing one term?

texts:
I am solving a question that says if we have sequence $x(n)$ of a signal missing one term then we have to find autocorrelation matrix $R_x$ as follows:

$$R_x = E\{\mathbf {xx^H}\}$$
Now if I take $x(n) = [x(0) , x(2),x(3)]$ that is missing term $x(1)$
I found following value of $R_x$

$R_x =  $$\left(
        \begin{matrix}
        r_x(0) & r_x^*(2) & r_x^*(3) \\
        r_x(2) & r_x(0) & r_x^*(1) \\
        r_x(3) & r_x(1) & r_x(0) \\
        \end{matrix} \right)
$

Now I check that this matrix is Toeplitz or not, by the definition of Toeplitz (constant diagonal term) $R_x$ is Toeplitz. But when I checked the solution it says the above matrix is not Toeplitz and there is no explanation given. I wonder what I am missing in definition of a Toeplitz matrix that does not fit in the above matrix. Can anyone explain this?


texts: How to do an FFT of a line profile obtained by imagej?

texts:
First of all, I'm a total "signal processing" newbie. So, bear with me.

I know that in ImageJ/Fiji I can use the Plot Profile function to select a line profile of an image. Is there a simple way to plot a 1D FFT transform of the data along the line?

If there is a better SO site for this kind of question, please let me know.


texts: Scipy stft, how is shape of f determined?

texts:
Using scipy stft, like this :

sample_rate, samples = scipy.io.wavfile.read(sample_dir/'original'/test_sample_name)
f, t, Zxx = scipy.signal.stft(samples.transpose(), sample_rate)

f.shape is equal to 129. My understanding is that this number represents the number of frequency at which we evaluate the Fourier Transform. How is it calculated? Can we fix this number?


texts: Drawing the modulus from a Transfer function

texts:
As you may guess from my other questions i am a student so pardon me any ignorance and guide me to the truth. Thanks

In this exercise from an exam i was given that transfer function with poles and zeros. I was asked to draw the modulus and the "UHD imagive below the Transfer function" image was my answer.

I know that when you encounter a zero, you "go down" and when you encounter a pole you "go up" in the graph. When you encounter both "aligned", you get a cero as they cancel each other.

The answer from the teacher was that the maximum, which is 1, is at pi/8 (and -pi/8), not at pi/4, because at pi/4 (and -pi/4) he draw the minimum.

Why is that? and what i did got wrong or missed in my understanding?

Thanks for answer.

EDIT: Here i add another example from a recent exam. I would like if someone could do a brief/short explanation in how to properly draw it and what should i know or notice.

The question is the same, but i lack proper understanding of it.


texts:
I have a set of $n$ low frequency time-varying signals $y(t)$.

I would like to express the variability within this set of signals using a single number. My current idea is to:

Compute the average signal $\bar y(t)$ for the set using

$$\bar y(t) = \dfrac{y_1(t) + y_2(t) + \ldots + y_n(t)}{n}  $$
Compute the 'standard deviation' of each individual signal with respect to the mean signal $\bar y$ using the Pearson correlation $r$

$$\bar r = \dfrac{ r(y_1(t), \bar y(t)) + r(y_2(t), \bar y(t)) + \ldots + r(y_n(t), \bar y(t))}{n}$$

This would then give me a sense of the variability within the set of signals.

My question is whether there are better, different, or more established ways of doing this?


texts: Stopband's edges in bandpass IIR filter in MATLAB

texts:
I've tried to get an IIR bandpass filter having the passband [0.1,0.5] Hz and working at the sampling frequency 8 Hz using the following command:
[fdata, BPDigFilt] = bandpass(data, [0.1,0.5], 8,        ...
                             'ImpulseResponse', 'iir',   ...
                             'StopbandAttenuation', 30);

I get the following frequency response:

but the command info(BPDigFilt) returns this
Second Passband Edge   : 500 mHz
Passband Ripple        : 0.1 dB
First Stopband Edge    : 84.3 mHz
First Passband Edge    : 100 mHz
First Stopband Atten.  : 30 dB
Second Stopband Atten. : 30 dB
Second Stopband Edge   : 1.0495 Hz

The strange thing is the second stopband edge equal to 1.0495 Hz. It seems strange to me because, from the previous plot, it seems that the required attenuation is reached at 0.6 Hz about, and not near 1 Hz.
Can someone explain to me the reason of this behaviour?


texts: Why do I have frequency leakage in DFT after zero padding if frequency resolution is fine?

texts:
Let's consider this example:

Fs=1000;
Ns=500;
t=0:1/Fs:(Ns-1)*1/Fs;
f1=10;
f2=400;
x=5+5*sin(2*pi*f1*t)+2*sin(2*pi*f2*t);
X=fft(x);

In this scenario, frequency resolution is 2, and all frequency components are captured correctly.
However, if I do this:

  X=fft(x,1000);

frequency resolution is 1, but there is a spectral leakage. Similar effect is seen here.
It seems to me, that Fourier transforms of both windows (one with length 500, and one with length 1000) have zeros at the frequencies that are presented in the signal, so I don't see why leakage will happen?


texts:
I'm trying to design a digital low pass filter with a narrow transition band. My sampling rate is 25 kHz, the cut off frequency is 60 Hz & the transition band width is 4 Hz. I'm looking for about 40 dB attenuation in the stop band & 0.1 dB in the pass band.

I looked at the traditional FIR methodology using the Windowed Sinc filter but the number of taps required to get such a transition band is too high. Going by the reference here, we get a value of 25000 which is 1 second in terms of time.

As an alternative, I looked into IFIR filters. While the model filter results in a decrease in the filter length by a value that's equal to the interpolation factor, the interpolated filter formed by inserting zeros into the model filter's kernel increases the length to almost the same as what you'd get with a conventional FIR filter.

So, my question is, does the IFIR methodology make the filter more efficient from a group delay point of view? Is the only benefit here the fact that more zeros implies fewer multiplications & additions? I'm trying to achieve around 1 to 2 cycles of delay (in terms of time for a 50 Hz signal). Considering that the sampling rate is so high, shouldn't it be easier to achieve this? If not, what other methods could be used to obtain such a filter?


texts: Why does an anti-symmetric function has zero amplitude at the center of an even length window

texts:
I am performing FFT on a real odd function and the resultant transform has zero amplitude in the last bin. Essentially if Y= rfft(X), then Y[-1] is always zero. I stumbled on this answer which says

  An anti-symmetric waveform has to be zero at the center of an even length window.

Now, this might be very easy to prove but I cannot think of why this is true. If I use the formula then to the first term should be zero, not $N/2$ th. It's clear to me why it should be real but why should the amplitude be zero for real-valued, odd function. Any help will be appreciated.


texts: How much zero padding I should do to an audio signal before fft?

texts:
I was working on project and I need to do fft to my audio signal. I was going through a code and found following line. Can anyone explain me the line of code.
padlen = int((numframes - 1) * frame_step + frame_len)
zeros = numpy.zeros((padlen - slen,))

Can anyone explain  above line? (padlen)
Here is the complete code
slen = len(sig)
    frame_len = int(round_half_up(frame_len))
    frame_step = int(round_half_up(frame_step))
    if slen <= frame_len:
        numframes = 1
    else:
        numframes = 1 + int(math.ceil((1.0 * slen - frame_len) / frame_step))

    padlen = int((numframes - 1) * frame_step + frame_len)

    zeros = numpy.zeros((padlen - slen,))
    padsignal = numpy.concatenate((sig, zeros))

where
frame_len = sample_rate * frame_size(0.025s)
frame_step = sample_rate * step_size( 0.010s)
Other Code I found is
def frame_audio(audio, FFT_size=2048, hop_size=10, sample_rate=44100):
    # hop_size in ms

    audio = np.pad(audio, int(FFT_size / 2), mode='reflect')
    frame_len = np.round(sample_rate * hop_size / 1000).astype(int)
    frame_num = int((len(audio) - FFT_size) / frame_len) + 1
    frames = np.zeros((frame_num,FFT_size))

    for n in range(frame_num):
        frames[n] = audio[n*frame_len:n*frame_len+FFT_size]

    return frames

Any explanation about this piece.


texts: How to create a synthetic time series where power spectral density estimation is achieves better results than a direct Fourier transform?

texts:
I am trying to create a synthetic time series where PSD estimation is necessary and useful to recover the correct spectral information of the time series. But so far I can only create a time series where either a) Fourier transform is superior and recovers sharper peaks; or b) Neither Fourier transform or PSD recover useful spectral information.

I have been trying to show this by creating a time series in MATLAB which has a peak frequency of 40 Hz but is also a wide-sense stationary random process. I then want to use power spectral density estimation to recover the spectral content (with a peak at 40 Hz).

I have tried doing a variety of things including:

1) Create a time series and add random noise. Both the Fourier transform and the PSD recover the true spectrum and the Fourier transform recovers sharper peak. So why use PSD?

2) Create a time series and add a random phase. Fourier transform recovers sharper peak at 40 Hz. Why use PSD?

3) Create a time series with random noise and then randomly sample the time series by pulling data samples from the time series with replacement. This results in random spectra with no peaks using either Fourier transform or PSD.

4) Create a time series with random noise and then re-sample the time series at a lower sample rate. This results in both Fourier transform and PSD recovering the peak at 40 Hz but, once again, the Fourier transform has sharper peaks.

What am I doing wrong? How can I create a synthetic wide-sense stationary time series in MATLAB?


texts: The Effect of the Standard Deviation ($ \sigma $) of a Gaussian Kernel when Smoothing a Gradients Image

texts:
I am trying to understand why smoothing an image with a gaussian kernel of different sigma values and then computing the gradients of the smoothed image leads to "thicker" trails. In the image below, the left image is produced by convolving the image with the derivative of a gaussian kernel where $\sigma = 1$ while the right picture shows image gradients when the image is convolved with a gaussian kernel where a $\sigma = 2$.
From what I understand, the gaussian kernel with a higher $\sigma$ value tends to make pixels look "more" like it's neighbors, because it places higher weights on pixels that are further away from the middle.


texts: What kind of information can we obtain from curvature for exemplar-based inpainting?

texts:
I just want to know what kind of information does Curvature suggest to us? I've seen an article that used from the curvature of an Image. I can't understand what's the purpose of the author or it's better to say that what kind of information Curvature has!!!

Can anyone help me to handle this problem?

The link of the article:
Article link

The article defines isophote curvature as:
$$
K(p) = \bigtriangledown \left[ \frac{ \bigtriangledown I_p}{ \left|  \bigtriangledown I_p \right|} \right]
$$


texts: I am working with ultrasonic transducers with frequencies 100 kHz, one for transmission one for reception. The transmitted signals are chirps. I would like to measure the time of flight between transmission and reception. Can anyone help me with how could I do the phase/amplitude equalization of impulse response of transducer? My received signal is distorted and attenuated.


texts:
I am reading a paper on frequency warping and I need to do a little manipulation of the Z-transform. Can somebody help me on how can I go about deriving equations $(3)$ and $(4)$ from equations $(1)$ and $(2)$.

$$ \tag{1}
\bar{A}(z) = \sum_{k=0}^{p}\bar{a}_k\bar{z}^{-k}(z), \ \ \ \ \ \ \bar{a}_0=1
$$

$$ \tag{2}
\bar{z}^{-1} = \dfrac{(1-a^2)z^{-1}}{1-az^{-1}}-a
$$

putting $(2)$ in $(1)$ and using binomial formulas:

$$ \tag{3}
\bar{A}(z) = \sum_{k=0}^{p}b_k\tilde{z}^{-k}(1-az^{-1})^{-k}
$$

$$ \tag{4}
b_k=\sum_{n=k}^{p}C_{kn}\bar{a}_n, \ \ \ \ \ \ C_{kn}= {n \choose k}(1-a^2)^k(-a)^{n-k}
$$

Thanks in advance!